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ABSTRACT 

The  literature  has  shown  that  crack  propagation  in  cracked 
metal  sheets  can  be  significantly  reduced  by  bonding  an  uncracked 
reinforcement  to  the  metal  sheet.  However,  cyclic  debonding  typically 
occurs  over  a  localized  area  near  the  crack.  Herein,  an  analysis  was 
developed  to  predict  both  the  crack  growth  and  debond  growth  in  a 
reinforced  system.  The  analysis  was  based  on  the  use  of  complex 
variable  Green's  functions  for  cracked,  isotropic  sheets  and  uncracked, 
orthotropic  sheets  to  calculate  inplane  and  interlaminar  stresses, 
stress  intensities  and  strain-energy-release  rates.  An  iterative 
solution  was  developed  that  used  the  stress  intensities  and  strain- 
energy-  release  rates  to  predict  crack  and  debond  growths,  respectively, 
on  a  cycle-by-cycle  basis.  The  analysis  was  verified  with  experiments. 

The  analysis  was  used  in  a  parametric  study  of  the  effects  of 
boron-epoxy  composite  reinforcement  on  crack  propagation  in  aluminum 
sheets.  The  study  showed  that  the  size  of  the  debond  area  has  a 
significant  effect  on  the  crack  propagation  in  the  aluminum.  For 
small  debond  areas  the  crack  propagation  rate  is  reduced  significantly. 
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but  these  small  debonds  have  a  strong  tendency  to  enlarge.  Debond 
growth  is  most  likely  to  occur  in  reinforced  systems  that  have  a 
cracked  metal  sheet  reinforced  with  a  relatively  thin  composite  sheet. 

The  analysis  predicts  crack  growth  in  reinforced  systems. 

Hence,  the  analysis  can  be  applied  in  developing  methods  to  repair 
damaged  metal  structures  and  to  increase  the  lives  and  payloads  of  metal 
structures  by  selective  reinforcement. 
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CHAPTER  I 


THE  USE  OF  COMPOSITE  REINFORCEMENT 
TO  PREVENT  FATIGUE  FAILURE  IN  AIRCRAFT 

A  potential  cause  of  aircraft  crashes  is  fatigue  failure.  As 
shown  by  Hardrath  (1971),  most  types  of  civil  aircraft  have  experienced 
some  form  of  fatigue  problem.  In  addition  Lowndes  and  Miller  (1969) 
indicate  that  fatigue  failures  have  frequently  occurred  in  military 
aircraft.  In  some  cases  the  fatigue  failures  led  to  loss  of  lives  and 
the  aircraft.  In  efforts  to  eliminate  such  failures,  both  government 
research  laboratories  and  aircraft  manufacturers  have  studied  the 
fatigue  failure  process  in  depth.  These  studies  showed  that  the  rate 
at  which  the  fatigue  damage  develops  in  metals  is  a  function  of  the 
stress  level  in  the  structure  and  occurs  in  three  stages:  crack 
initiation,  stable  crack  propagation,  and  unstable  crack  propagation 
(catastrophic  failure).  Although  aircraft  structures  can  be  designed 
to  have  stresses  low  enough  to  prevent  fatigue  failures,  the  weight 
penalty  would  be  enormous  and  would  make  the  aircraft  uneconomical 
to  operate.  Hence,  a  trade-off  exists  between  low  stresses  and  low 
weight,  and  weight  efficient  structures  will  almost  always  have 
stresses  high  enough  to  support  fatigue  damage  accumulation. 

Fatigue  cracks  initiate  at  local  stress  concentrations  in  the 
structure.  The  local  stress  concentrations  may  be  caused  by  poor 
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fatigue  design,  by  manufacturing  defects  in  the  material,  or  by  damage 
caused  by  the  flight  environment.  Although  methods  can  be  employed  to 
reduce  the  occurrence  of  fatigue  crack  initiation,  the  development  of 
such  cracks  seems  almost  inevitable. 

Once  a  crack  initiates  it  grows  at  a  stable  rate  until  it  reaches 
some  predictable,  critical  length  after  which  catastrophic  failure  follows 
Fortunately,  in  aluminum  aircraft  structures  the  critical  crack  length 
is  large  and  the  crack  is  easy  to  detect  long  before  it  reaches  a 
critical  length.  Consequently,  fatigue  cracks  can  be  tolerated  in 
an  aircraft  structure  as  long  as  the  structure  is  inspected  periodically 
to  locate  cracks  before  they  become  critical.  Of  course,  once  the 
crack  is  detected  it  must  be  repaired  before  it  becomes  critical.  The 
repair  can  be  made  by  either  replacing  the  component  or  by  repairing 
it  in  situ.  Because  replacing  a  component  may  involve  high  cost  and 
keep  the  aircraft  out  of  service  for  extended  periods  of  time,  repair¬ 
ing  the  component  in  situ  is  frequently  very  desirable. 

Basically,  a  fatigue  crack  can  be  repaired  by  reducing  the  stress 
state  in  the  vicinity  of  the  crack  tip.  One  method  of  reducing  the 
stress  state  is  to  reinforce  the  crack  with  unidirectional  composite 
(fibers  are  perpendicular  to  plane  of  crack).  The  composite  reinforce¬ 
ment  reduces  the  stress  state  near  the  crack  tip  by  two  mechanisms. 

First,  adding  the  composite  reinforcement  lowers  the  overall  stress  in 
the  cracked  metal  by  increasing  the  cross-sectional  area  and  by 
providing  an  alternate,  stiffen  load  path  by  virtue  of  the  high 
modulus  of  the  composite.  This  reduction  in  stress  can  be  easily 
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calculated  by  simple  strength  of  materials  theory  and,  hence,  is 
easy  to  investigate.  The  second  mechanism  conies  from  the  development 
of  stresses  between  the  metal  and  composite  adherends.  Several  papers 
(Kula  et  al  1973,  Ellis  1976,  Johnston  and  Stratton  1975,  Ratwani  1977) 
have  shown  that  these  interlaminar  stresses  have  a  profound  effect  on 
crack  propagation.  These  interlaminar  stresses  reduce  the  stresses 
near  the  crack  tip  and,  consequently,  retard  its  growth.  As  will  be 
shown  later,  the  investigation  of  these  interlaminar  stresses  requires 
a  much  more  extensive  analysis  than  that  provided  by  strength  of 
materials  theory. 

Consequently,  a  need  exists  for  the  development  of  a  realistic 
fatigue  analysis  that  incorporates  the  effects  of  the  interlaminar 
stresses.  Accordingly,  the  objective  of  this  dissertation  was  to 
develop  such  an  analysis  and  use  it  to  study  the  effects  of  composite 
reinforcement  on  the  fatigue  life  of  cracked  metallic  structure. 

To  meet  this  objective  the  following  approach  was  taken.  First, 
in  Chapter  II  the  fatigue  behavior  of  the  constituents  of  the  reinforced 
system  was  characterized.  Next  in  Chapter  III  the  fatigue  behavior  of 
the  reinforced  system  was  studied  experimentally.  Then,  in  Chapter  IV 
with  the  use  of  the  results  of  Chapters  II  and  III  and  complex  variable 
theory,  a  static  analysis  was  developed  that  related  applied  loads, 
adherend  thicknesses,  debond  size,  and  crack  length  to  crack  propagation 
rates.  Next,  in  Chapter  V  the  analysis  was  further  developed  to  predict 
both  debond  and  crack  growth  as  a  function  of  applied  load  cycles.  The 
accuracy  of  the  analysis  was  investigated  in  Chapter  VI.  Finally,  in 
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Chapter  VII  the  analysis  was  used  to  parametrically  study  crack 
propagation  in  reinforced  systems. 

Throughout  the  development  of  the  analysis  many  items  required 
detailed  analytical  or  experimental  investigation.  These  investigations 
were  developed  and  discussed  in  several  appendices. 


CHAPTER  II 


FATIGUE  BEHAVIOR  OF  THE  REINFORCED  SYSTEM  CONSTITUENTS 

The  reinforced  system  that  will  be  considered  herein  is  composed 
of  adherends  made  out  of  two  dissimilar  materials,  an  aluminum  sheet 
and  a  composite  sheet,  that  are  bonded  to  each  other  with  a  relatively 
thin,  room- temperature  curing  adhesive.  This  system  is  intended  to 
model  the  repair  of  a  cracked,  aluminum  aircraft  component  that  is 
repaired  by  bonding  a  composite  sheet  to  it.  Each  of  the  three 
constituents  of  the  system  -  the  metal,  the  composite,  and  the  adhesive 
exhibits  different  fatigue  behavior  and  plays  an  important  role  in  the 
fatigue  behavior  of  the  reinforced  system.  Consequently,  to  analyze 
the  fatigue  process  in  the  system,  the  fatigue  behavior  of  each  of  the 
constituents  needs  to  be  understood.  In  the  following  sections  the 
fatigue  behavior  of  each  of  the  constituents  will  be  discussed  and 
analysis  methods  formulated. 

Fatigue  of  Metals 

As  pointed  out  by  Erodogan  (1968),  the  fatigue  process  in  metals 
occurs  in  three  different  stages:  crack  initiation,  stable  crack 
growth,  and  unstable  crack  growth  (fracture).  Current  aircraft  design 
methods  focus  on  the  latter  two  stages  of  the  fatigue  process  by  using 
a  "Damage  Tolerant  Design  Philosophy"  (military  specification 
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MIL-A-83444) .  This  philosophy,  as  far  as  fatigue  damage  accumulation  is 
concerned,  admits  that  initial  flaws  such  as  cracks,  exist  in  aircraft 
components  that  are  fatigue  critical,  i.e.  may  fail  under  cyclic  loading. 
But,  the  philosophy  also  assumes  that  these  initial  cracks  grow  stably 
and  can  be  detected  during  periodic  inspections  before  they  reach  a 
critical  crack  length.  Once  the  damage  is  detected,  it  can  be  repaired. 
Hence,  the  validity  of  this  philosophy  rests  on  the  accurate  prediction 
of  crack  growth  rate  and  critical  crack  length.  Fracture  mechanics 
theory  can  be  used  to  determine  both  crack  growth  rate  and  critical 
crack  lengths. 

Fracture  mechanics  theory  was  conceived  when  Griffith  (1921) 
related  fracture  to  an  energy  balance  as  the  crack  extended.  In  1957 
Irwin  related  the  stress  state  at  the  crack  tip  to  fracture.  A 
schematic  of  a  crack  tip  and  equations  for  the  stresses  very  close  to 
it  (Sih  and  Liebowitz  1968)  are  shown  on  figure  1  (additional  terms  not 
shown  in  the  equations  have  a  negligible  effect  on  the  stress  state 
near  the  crack  tip).  As  may  be  seen  from  the  equations,  as  the  distance 
from  the  crack  tip,  r,  approaches  zero  the  stresses  become  infinite. 
Consequently,  at  the  crack  tip  where  the  stresses  are  infinite  a 
singularity  exists.''  The  coefficients  of  the  stress  distributions, 
ki  and  are  the  stress  intensity  factors  which  are  used 

extensively  in  fracture  mechanics.  The  Mode  I  stress  intensity 
designated  by  kj  is  associated  with  the  stresses  that  deform 
the  crack  surfaces  symmetrically  with  respect  to  the  original  plane  of 

^In  reality  infinite  stresses  cannot  exist  in  the  material  and 
local  yielding  of  the  material  occurs.  This  local  yielding  is  ignored 
in  linear  elastic  fracture  mechanics. 


8 


the  crack  while  the  Mode  II  stress  intensity  designated  by  ka  is 
associated  with  stresses  that  cause  shear  displacements  between  the 
crack  surfaces.  In  1960  Sanders  showed  that  the  stress  intensities 
were  related  to  the  strain  energy  release  as  the  crack  extended.  Hence, 
the  stress  state  near  the  crack  tip  was  related  to  the  Griffith  theory 
of  fracture,  and  the  foundation  of  fracture  mechanics  was  formed. 

The  stress  intensities  can  be  determined  from  complex  stress 
functions  determined  from  the  theory  of  elasticity  (Sih  and  Liebowitz 
1968)  as 

ki  -  ik2  =  2/2  lim  {/z^I  (z)}  (1) 

z  a 


where 


z  =  X  +  iy  and  i  =  /T 

and  X  and  y  are  the  cartesian  coordinates  and  $(z)  is  the  complex 
stress  function  as  developed  by  Muskhelishvili  (1975)  that  satisfies 
the  equations  (plane  strain  or  stress) 

a  +  0  =  2[$(z)  +  $(z) ]  (2) 

X  y 


^i^xy  ■  ^x  '^y  ^  2[z<^'(z)  +  'V{z)]  (3) 

where  a  ,a  and  a  are  the  stresses  in  the  cartesian  coordinates 
X  y  xy 

and  H'(z)  is  another  stress  function.  The  two  stress  functions  are 
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functions  of  both  the  loading  conditions  and  the  configuration  of  the 

body  and  will  be  discussed  in  detail  in  later  chapters. 

The  stress  intensities  can  be  related  to  both  crack  propagation 

rates  and  critical  crack  lengths.  On  the  basis  of  the  Griffith  theory 

of  fracture,  a  critical  value  of  the  strain  energy  release  can  be  found 

and  hence,  according  to  Sanders  (1960)  a  critical  value  of  the  stress 

intensity  can  be  found.  For  the  material  used  in  this  study,  2924-T3 

aluminum,  Hudson  (1969)  showed  that  the  critical  value  for  k  is 

ic 

56,000  psi-in  .  Hence,  with  the  use  of  equation  (1)  and  the  appropriate 
stress  functions,  the  fracture  can  be  predicted. 

Cyclic  crack  growth  rates  were  related  to  the  stress  intensity 
by  Paris  (1961)  by  the  empirical  formula 

da/dN  =  C(k  )** 

1 

where  da/dN  is  the  crack  propagation  rate,  C  is  am  empirical 
constant  and  k^  indicates  the  stress  intensity  range  during  cyclic 
loading.  Forman  et  al  (1967)  improved  this  equation  by  including  the 
critical  stress  intensity  k^^  and  the  stress  ratio  R,  which  is 
the  ratio  of  the  minimum  to  maximum  stress  in  the  loading  cycle, 
in  the  empirical  formula 

c^(k^)"i 

da/dN  =.  - ^ ^ -  (4) 

(l-R)k  .  -  k 

IC  1 

where  c  and  n  are  empirical  constants  and  =  56,000  psi-in"^. 
11 
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For  2024-T3  aluminum  Hudson  (1969)  showed  that  the  constants  and 
n^  have  the  values 

=  3.22  X  10'^^ 

n^  =  3.38 

As  the  previous  discussion  implies,  once  the  stress  functions  for 
a  cracked  body  are  known,  the  stress  intensities  can  be  calculated. 

With  the  stress  intensities  both  the  crack  propagation  rate  and  critical 
crack  length  can  be  predicted.  The  crack  propagation  rate  and  critical 
crack  length  can  be  used  in  a  Damage  Tolerant  Design  Philosophy  to 
predict  life  of  aircraft  components.  The  life  is  predicted  by  first 
assuming  that  the  structure  contains  cracks.  The  lengths  of  these 
cracks  are  assumed  to  be  the  largest  crack  detected  in  the  structure 
or  the  largest  crack  which  can  be  overlooked  due  to  the  resolution 
of  the  inspection  technique.  Then,  by  using  the  assumed  or  detected 
crack  length  and  fracture  mechanics  theory,  the  number  of  load  cycles 
to  fracture  can  be  predicted.  On  the  basis  of  these  calculations, 
inspection  intervals  are  determined  to  assure  that  cracks  can  be 
detected  and  repaired  before  they  reach  a  critical  length. 

Fatigue  of  Bonded  Systems 

To  perform  a  realistic  fatigue  analysis  of  the  reinforced  system, 
the  fatigue  behavior  of  the  adhesive  in  situ,  herein  called  the  bond, 
must  be  characterized.  Several  researchers  have  shown  that  the  bond 
deteriorates  when  subjected  to  a  cyclic  load.  Within  this  disserta¬ 
tion  this  deterioration  will  be  called  debonding.  Hoffman  and  June 
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(1973)  studied  debonding  by  recording  the  debond  propagation  as  a  func¬ 
tion  of  applied  load  cycles.  They  showed  that  a  myriad  of  factors 
such  as  type  of  adhesive,  adherend  and  adhesive  thickness,  method  of 
curing,  and  aging  all  affect  the  debonding.  Roderick  et  al  (1976) 
showed  that  debonding  could  occur  as  failure  within  the  adhesive  as  a 
cohesive  failure,  between  the  adhesive  and  an  adherend  as  an  adhesive 
failure,  or  in  the  composite  material.  Because  of  the  variety  of 
failure  modes,  the  analysis  of  the  debonding  is  difficult.  The  most 
progress  in  analysis  of  debonding  appears  to  stem  from  the  energy 
approach  developed  by  Griffith  (1921). 

The  first  application  of  the  energy  approach  appears  to  be  by 
Rippling  et  al  (1964)  in  the  study  of  fracture  toughness  of  bonded 
joints.  Since  Rippling's  paper,  Mostovoy  and  Ripling  have  published 
several  other  papers  on  fracture  toughness  of  bonds:  Mostovoy  and 
Rippling  1966,  Mostovoy  et  al  1967,  Mostovoy  and  Rippling  1971.  However, 
a  correlation  between  the  fracture  energy  and  the  stress  state  near 
the  debond  tip  has  not  been  made  in  the  bonded  systems.  Wang  et  al 
(1976)  showed  that  a  primary  reason  for  the  lack  of  correlation  appears 
to  be  the  development  of  large  regions  of  plastic  yielding  in  the 
adhesive.  Hence,  linear  elastic  fracture  mechanics  based  on  small 
yield  zones  and  stress  intensities  at  a  crack  tip  do  not  appear 
applicable  to  bonded  systems. 

However,  by  applying  an  energy  approach,  Roderick  et  al  (1975) 
showed  that  the  debond  propagation  rates  can  be  correlated  for  specimens 
with  different  thickness  adherends  with  a  Paris  (1961)  type  equation 
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db/dN  =  C2(G)'^2  (5) 

where  both  and  ai^e  empirical  constants  for  a  specific  bond 
system  and  G  is  the  strain  energy  released  as  the  debond  extends. 

As  shown  by  Roderick  et  al  (1976),  the  parameters  c^  and 
x\2  vary  for  different  bonded  systems.  The  bond  system  used  in  this 
dissertation  was  2024-T3  aluminum  bonded  to  unidirectional  boron/epoxy 
with  a  room  temperature  curing  adhesive.  Shell  EA-934.  For  this  system 
the  empirical  constants  were  determined  by  methods  discussed  in 
Appendix  A  and  were  found  as 

=  3.158  X  10'^ 

nj  =  3.616 

With  these  constants,  a  value  for  the  strain  energy  release  rate,  G, 
and  equation  (5),  the  debond  growth  rate  can  be  predicted.  The 
calculation  of  G  for  debonding  in  the  reinforced  system  will  be 
discussed  in  detail  in  Chapter  V. 

Fatigue  of  Composite  Materials 

The  term  "composite"  may  refer  to  a  myriad  of  systems  composed 
of  a  wide  spectrum  of  different  types  of  fibers  and  matrices.  Further¬ 
more,  each  system  may  have  widely  different  fatigue  characteristics 
depending  upon  the  fiber  orientations,  stacking  sequences,  and  loading 
conditions.  Durchlaub  and  Freeman  (1974)  showed  that  fatigue  damage 
in  composites  may  occur  perpendicular  to,  parallel  to,  or  at  an 
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angle  to  the  loading  axis  depending  upon  the  fiber  orientation.  Foye 
and  Baker  (1970)  showed  that  the  lives  of  composite  laminates  could 
vary  as  much  as  an  order  of  magnitude  by  changing  their  stacking 
sequences.  Reifsnider  et  al  (1974)  showed  that  changing  the  frequency 
of  the  applied  cyclic  load  affects  both  the  mode  of  failure  and  the 
fatigue  life.  As  evident  from  these  observations,  the  fatigue  behavior 
of  composite  materials  is  complex. 

Currently  the  understanding  of  the  fatigue  process  in  composites 
appears  primitive  although  some  progress  in  developing  an  understanding 
has  been  made.  As  pointed  out  by  Salkind  (1973),  fatigue  failure  in 
composites  can  occur  in  different  failure  modes  such  as  matrix  cracking, 
delamination,  and  fiber  fracture.  Also,  evidence  exists  that  suggests 
that  the  fatigue  process  is  a  result  of  primarily  matrix  deterioration 
(Roderick  and  Whitcomb  1977).  If  matrix  deterioration  is  the  primary 
cause  of  fatigue  in  composites,  then  the  various  failure  modes  could 
be  explained  in  terms  of  different  stress  states  in  the  matrix  depend¬ 
ing  upon  the  fiber  orientation  and  stacking  sequences  of  a  specific 
laminate.  Hence,  those  laminates  in  which  the  matrix  is  highly  stressed 
would  most  likely  degrade  under  cyclic  load  while  those  laminates  in 
which  the  matrix  is  lightly  stressed  would  not. 

Following  this  line  of  thought,  composites  in  which  the  fibers 
transmit  the  load,  fiber  controlled  composites,  would  have  long  fatigue 
lives  while  those  in  which  the  matrix  transmits  the  load,  matrix 
controlled  composites,  would  have  short  lives.  An  example  of  a  fiber 
controlled  composite  is  a  unidirectional  laminate  loaded  along  the  axis. 
On  the  other  hand,  an  example  of  a  matrix  controlled  laminate  is  one 
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in  which  the  fibers  are  at  45°  angles  to  the  loading  axis.  As  shown  by 
Durchlaub  and  Freeman  (1974),  the  matrix  controlled  laminate  does  degrade 
rapidly  under  cyclic  loading  while  the  fiber  controlled  laminate  does  not. 

Two  basic  approaches  to  predict  the  diverse  fatigue  behavior  of 
composite  materials  are  currently  being  developed  by  researchers.  In 
the  first  approach,  the  laminate  behavior  is  described  by  a  statistical 
model  (Halpin  et  al  1972)  that  relates  the  static  strength  and  fatigue 
life  distributions  by  assuming  that  the  residual  strength  of  the 
laminates  degrades  monotonically  (Yang  and  Liu  1977).  Because  this 
method  is  based  on  experimental  results,  it  can  be  applied  easily. 
However,  it  does  not  apply  to  laminates  whose  residual  strength  does 
not  decrease  monotonically^  with  applied  load  cycles.  Also,  this 
method  requires  extensive  testing  every  time  the  stacking  sequence  or 
fiber  orientation  changes. 

The  second  approach  as  developed  by  McLaughlin  et  al  (1975) 
couples  basic  fatigue  data  on  the  laminae  level  with  a  stress  analysis 
to  predict  both  the  mode  of  fatigue  failure  and  the  fatigue  lives  of 
laminates.  Because  this  approach  is  based  on  laminae  data  rather  than 
laminate  data,  it  can  be  used  to  predict  the  fatigue  behavior  of  lami¬ 
nates  with  different  stacking  sequences  and  fiber  orientations  without 
extensive  testing.^  The  major  drawback  to  this  approach  is  its 

^Durchlaub  and  Freeman  (1974)  showed  that  the  residual  strength  of 
notched  laminates  could  increase  after  fatigue  loading. 

^The  analysis  originally  proposed  by  McLaughlin  et  al  did  not 
consider  interlaminar  stresses  and  therefore  could  not  account  for 
changes  in  stacking  sequences,  but  incorporation  of  interlaminar  stress¬ 
es  into  the  analysis  has  been  done  and  will  be  shown  in  a  NASA 
contractors  report  released  in  1978. 
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complexity  in  attempting  to  develop  simple,  realistic  stress  analyses 
and  a  failure  criterion  on  the  laminae  level. 

The  state  of  the  art  of  fatigue  analysis,  in  the  author's  opinion, 
is  still  in  the  early  stages  of  development  and  not  yet  capable  of 
reliable  life  predictions  for  general  laminates.  As  a  consequence 
the  fatigue  behavior  of  the  unidirectional  boron/epoxy  used  in  the 
present  study  cannot  be  described  by  relatively  simple  fatigue  analyses 
as  was  the  case  for  the  cracked  metal  sheet  and  the  bond  system.  Accord¬ 
ingly,  the  fatigue  behavior  will  be  determined  solely  by  experimental 
data. 

Shockey  et  a1  (1970)  showed  that  unidirectional  boron/epoxy 
laminates  that  were  loaded  along  the  fiber  axis  had  an  average  ultimate 
tensile  strength  of  193  ksi;  when  these  laminates  were  cycled  under 
constant  amplitude  cyclic  loading  with  R  =  0.1,  they  retained 
73  percent  of  their  ultimate  tensile  strength  after  10^  applied  load 
cycles.  Consequently,  in  an  attempt  to  prevent  failure  in  the 
unidirectional  boron/epoxy,  stress  along  the  fiber  axis  (based  on 
laminate  analysis)  was  kept  below  140  ksi. 

Having  discussed  the  fatigue  behavior  of  the  constituents  of  the 
reinforced  system  in  this  chapter,  the  next  chapter  deals  with  the 
fatigue  behavior  of  the  constituents  in  situ  in  the  reinforced  system. 
Hence,  the  next  chapter  discusses  fatigue  tests  of  reinforced  systems. 


CHAPTER  III 


FATIGUE  TESTS  OF  REINFORCED  SYSTEMS 

To  determine  the  fatigue  behavior  of  the  reinforced  system,  two 
large  panels  were  manufactured  and  tested.  The  panels  shown  on 
figure  2  were  made  of  8  x  24  inch  sheets  of  2024-T3  aluminum  and 
unidirectional  boron/epoxy.  EA-934  room  temperature  curing  adhesive 
was  used  to  join  the  sheets  with  the  bonding  process  described  in 
Appendix  A.  The  primary  difference  betvieen  the  panels  labeled  A  and  B 
on  figure  2  was  the  thickness  of  the  metal  and  composite  adherends.  To 
simulate  a  crack,  the  metal  adherend  contained  a  through-the-thickness 
narrow  slit  0.01  inch  wide  and  2  inches  long.  The  slit,  which  was  made 
by  an  electrical  discharge  process,  was  centered  along  the  horizontal 
centerline  of  the  panels.  In  both  panels  the  fibers  of  the  unidirec¬ 
tional  composite  run  parallel  to  the  longitudinal  axis  of  the  panels. 

The  panels  were  tested  in  a  300,000  pound  load  capacity  servo- 
hydraulic  fatigue  machine.  Both  panels  were  tested  under  constant 
amplitude  loading  with  R,  the  ratio  of  the  minimum  to  maximum  stress 

4 

in  the  load  cycle,  equal  to  0.01  at  a  test  frequency  of  2.5  Hz. 

For  the  fatigue  tests  of  both  panels  the  distance  between  test 
machine  grips  was  16  inches.  The  maximum  loads  applied  to  each  panel 

^The  test  frequency  was  limited  to  2.5  Hz  instead  of  the  10  Hz 
used  to  characterize  the  debond  behavior  in  Appendix  A  because  of 
test  machine  limitations. 
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.  2.  Panel  configuration 
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during  the  fatigue  tests  and  the  corresponding  stresses  in  the  adherends 
calculated  from  membrane  laminate  theory  (see  Appendix  C)  are  shown 
below. 


Panel 

Maximum  Load 

Stress, 

psi 

lbs 

metal 

compos i te 

A 

37,600 

19,600 

58,600 

B 

22,500 

14,600 

43,100 

During  the  fatigue  tests,  crack  lengths  were  measured  periodically 
with  an  optical  microscope.  Table  1  shows  the  measured  crack  lengths 
and  applied  load  cycles  for  both  tests.  The  crack  lengths  are  plotted 
against  the  applied  load  cycles  on  figure  3.  Note  that  on  the  figure 
the  abscissa  is  logarithmic  and  the  ordinate  starts  at  the  initial  half¬ 
crack  length  of  a  =  1.0  inch.  The  crack  propagation  rates  for  the 
panels  are  the  slopes  of  the  curves  shown  on  figure  3.  These  rates 
are  tabulated  in  Table  1  and  plotted  against  the  half-crack  length  on 
figure  4.  As  evident  from  figure  4,  the  crack  growth  rate  is  about 
two  orders  of  magnitude  larger  in  Panel  A  than  in  Panel  B. 

The  crack  propagation  rates  in  these  panels  is  a  function  of 
debonding  between  the  adherends.  If  the  adherends  were  completely 
debonded  the  crack  propagation  rate  would  be  much  larger  than  if  no 
debonding  occurred.  To  investigate  the  effect  of  debond  size  on 
crack  propagation  rate,  the  test  panels  were  examined  with  an  ultrasonic 
C-scan  (details  of  the  C-scan  method  are  discussed  in  detail  by 
McMaster  1963)  after  the  half-crack  length  grew  to  1.0  inch.  Figure  5 
shows  the  C-scans  of  the  panels.  On  the  figure  the  dark  parts  of  the 
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TABLE  1 


CRACK  LENGTHS 

AND  CRACK  PROPAGATION 

RATES 

Panel  A 

Panel 

B 

half-crack 

cycles 

crack 

cycles 

crack 

length 

propagation 

propagation 

rate 

rate 

a,  in 

*N 

da/dN 

X  10-5 

*N 

da/dN 

X  10-5 

1.05 

44,800 

0.236 

1.10 

1,570 

65,940 

0.312 

1.15 

7.57 

81,975 

0.217 

1.20 

3,220 

105,000 

0.33 

1.25 

9.02 

120,130 

0.369 

1.30 

4,440 

133,675 

0.311 

1.35 

10.8 

149,740 

0.328 

1.40 

5,540 

164,980 

0.339 

1.45 

11.6 

179,695 

0.321 

1.50 

6,790 

195,250 

0.277 

1.55 

12.0 

213,310 

0.299 

1.60 

7,790 

230,025 

0.275 

1.65 

12.2 

248,180 

1.70 

8,690 

.... 

1.75 

13.0 

279,795 

0.295 

1.80 

9,690 

296,765 

0.323 

1.85 

13.0 

312,260 

0.339 

1.90 

10,680 

326,990 

0.312 

1.95 

342,990 

2.00 

•  .  •  • 

*N  -  Instead  of  listing  the  number  of  cycles  that  caused  crack 
growth  at  both  crack  tips,  the  average  number  of  cycles  is  given. 


N,  applied  load  cycles 
Crack  growth  versus  applied  load  cycles 


Panel  A 

o  o  ° 


O  O  O 


o°°°oOo  o 

Panel  B 


1.4  1.6  1.8 

half-crack  length 
a,  inches 


Experimental  crack  propagation  rates 
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grips 


Specimen  A 


Specimen  B 


F  =  -^  =  0.75 


F  =  ^  =  0.125 


Fig.  5.  Ultrasonic  C-scans  of  fatigued  specimens 
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C-scans  are  areas  where  debonding  has  occurred.  As  evident  from  the 
figure  debonding  occurred  over  an  elliptical  area.  Ratwani  (1977) 
has  observed  similar  elliptical  debonds  in  metal  laminates.  The  major 
axes  of  the  elliptical  debonds  were  nearly  equal  to  the  crack  length 
in  the  metal  adherend  of  the  panels.  As  measured  from  figure  5  the 
minor  axes  of  the  debonds  were  3.0  inches  and  0.50  inch  for  Panels  A 
and  B  respectively.  However,  when  the  C-scans  shown  on  figure  5  were 
made,  they  were  distorted  along  the  longitudinal  axes  of  the  panels.  By 
taking  into  account  this  distortion,  the  minor  axes  of  Panels  A  and 
B  v/ere  found  to  be  4.0  and  0.67  inches  respectively.  With  these 
corrected  values  of  debond  length  along  the  minor  axes,  the  ratios  of 
minor  to  major  axes  of  the  debonds,  F  =  b/a,  were  found  as  1.0  and 
0.14  for  Panels  A  and  B  respectively. 

The  experiments  discussed  in  this  chapter  showed  that  the  fatigue 
in  reinforced  systems  occurs  as  col  linear  crack  growth  and  debond 
growth  over  an  approximately  elliptical  area.  In  Chapters  IV  and  V 
an  analysis  will  be  developed  that  can  model  this  observed  behavior. 

In  Chapter  VI  the  analysis  will  be  verified  by  comparing  the  experimental 
results  of  this  chapter  with  results  of  the  analysis. 


CHAPTER  IV 


STRESS  ANALYSIS 

As  shown  in  th©  pr©vious  chapter,  under  cyclic  loads  the  rein¬ 
forced  system  exhibits  both  crack  and  debond  growth.  Intuitively, 
the  rate  at  which  the  debond  and  crack  propagates  is  a  function  of  the 
stress  state  and  level  in  the  system.  Consequently,  to  predict  these 
rates  a  realistic  stress  analysis  is  required.  For  the  stress  analysis 
to  be  realistic  it  must  predict  stresses  in  the  adhesive  as  well  as 
in  the  adherends  of  the  system.  Because  adhesives  typically  exhibit 
nonlinear  behavior  (Hughes  and  Rutherford  1968),  the  stress  analysis 
must  include  nonlinear  behavior  of  the  adhesive.  The  first  step  in 
development  of  a  realistic,  nonlinear,  elastic  stress  analysis  is  the 
formulation  of  a  linear  elastic  stress  analysis. 

Formulation  of  Linear  Elastic  Solution 

To  formulate  an  elastic  solution,  the  reinforced  system  shown  in 
figure  6  was  considered.  As  shown  in  figure  6,  the  system  consists  of 
a  cracked  metal  sheet  bonded  to  a  composite  sheet  with  an  elliptical 
debond  between  the  sheets.  The  system  was  subjected  to  a  remote 
stress,  s.  A  rigorous  stress  analysis  of  this  system  requires  a 
three-dimensional  formulation.  Although  a  general,  exact  solution  is 
not  available,  finite  element  or  finite  difference  numerical  methods 
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can  be  employed  to  obtain  a  tractable  solution.  However,  these  solutions 
are  not  efficient  for  analyzing  reinforced  systems  in  which  the  crack 
length  and  debond  area  are  continually  changing.  An  alternate,  simple 
analysis  can  be  developed  by  assuming  that  the  adherends  are  in  plane 
stress  while  the  adhesive  is  in  pure  shear.  These  assumptions  were 
first  used  in  an  analysis  by  Volkerson  in  1934. 

The  validity  of  these  assumptions  for  analysis  of  the  reinforced 
system  shown  in  figure  6  was  investigated  in  detail  in  Appendix  B. 

As  shown  in  Appendix  B  with  a  simple  example,  the  assumption  can  lead 
to  errors  as  much  as  100  percent  in  the  calculated  adhesive  stresses 
as  compared  to  more  rigorous  finite  element  solutions.  Evidently, 
significant  shearing  deformation  occurs  in  the  adherends  of  the 
reinforced  system.  The  presence  of  the  adherend  shearing  deformation 
violates  Volkerson 's  assumptions,  but  as  shown  in  Appendix  B  an 
effective  shear  modulus,  6^^^,  can  be  determined  and  used  with  the 
assumptions  to  calculate  adhesive  shearing  stresses  within  a  few  percent 
of  the  finite  element  results. 

Arin  and  Erodogan  (1972)  used  Volkerson's  assumptions  with  complex 
variable  elasticity  theory  developed  by  Muskhelishvili  (1953)  and 
Lekhnitski  (1956)  to  analyze  a  system  similar  to  that  shown  in  figure  6. 
The  linear  elastic  stress  analysis  developed  herein  basically  follows 
the  concepts  used  by  Arin  and  Erodogan,  but  differs  in  the  formulation 
of  the  Green's  functions  used  in  the  elasticity  solution,  the  method  of 
numerical  integration  of  the  Green's  functions,  and  the  domain  of 
integration.  To  develop  the  stress  analysis,  the  reinforced  system 
is  free  bodied  as  shown  in  figure  7  (adhesive  layer  not  shown) 


Inplane  and  interlaminar  stresses  in  the  reinforced  system 
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using  Volkerson's  assumptions.  In  the  figure  the  remote  stress,  s, 
refers  to  the  applied  load  over  the  total  cross  sectional  area  of  the 
reinforced  system.  By  the  use  of  laminate  theory  as  described  in 
Appendix  C,  the  inplane  stresses  and  (where 

the  first  subscript  refers  to  stress  direction  and  the  second  subscript 
refers  to  the  metal  or  composite  adherend)  can  be  easily  calculated.  On 
figure  7,  f^^  f^^  indicate  shear  stresses  in  the  adhesive 

layer.  Throughout  this  dissertation  the  adhesive  shear  stresses  which 
will  be  assumed  to  act  as  body  forces  on  the  adherends  will  be  frequently 
called  interlaminar  stresses.  To  form  a  governing  equation,  these 
interlaminar  stresses  were  related  to  the  displacement  of  the  adherends 
in  the  following  manner. 

First,  the  shear  strain  in  the  adhesive  layer  was  related  to  the 
displacement  in  the  adherends  by  the  relations 


Y 


xz 


^m  -  ^^c 


ad 


(6) 


where  t  j  is  the  thickness  of  the  adhesive,  u  and  v  refer  to 
dd 

displacements  in  the  x  and  y  direction  respectively,  and  the 
subscripts  m  and  c  refer  to  the  metal  and  composite  adherends 
respectively.  Next,  Hook's  law  was  used  to  relate  the  shearing  strain 
in  the  adhesive  to  the  interlaminar  stresses  as 


T  =  GY 


(7) 
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Then,  by  the  use  of  the  effective  shear  modulus  and  equations  (6)  and 
(7)  the  interlaminar  stresses  can  be  related  to  the  adherend  dis¬ 
placements  as 


Equation  (8)  can  be  written  for  several  different  points  in  the 
reinforced  system  to  form  a  system  of  simultaneous  equations.  These 
simultaneous  equations  are  developed  in  detail  in  subsequent  sections 
of  this  chapter. 

The  displacements,  u  and  v,  in  the  adherends  were  related 
to  the  inplane  adherend  stresses  and  the  interlaminar  stresses  by 
several  functions  -  Fg  as 


The  displacements,  F^  and  Fj,  in  the  metal  adherend  due 
to  the  remote  stress  were  calculated  in  terms  of  two  stress  functions 
(!)(z)  and  w{z)  (Muskhelishvili  1975)  by  the  equation 
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2G(u^  +  iV|^)  =  ri<l)(z)  -  (jl)(z)  -  (z  -  z)$(z)  (10) 

3  -  V  _ 

where  n  = -  plane  stress,  i  =  /-I ,  z  =  x  +  iy,  the  bar 

1  +  V 

over  a  function  or  variable  denotes  the  complex  conjugate,  v 
is  Poisson's  ratio  and 


(13) 


In  equations  (11)  through  (13)  "a"  denotes  half  the  crack  length 

in  the  metal  sheet.  VJith  the  use  of  equation  (10),  the  functions 

and  F  were  found  to  be 
2 
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F  (am  ,am  )  =  Real 
1  XX’  yy' 


n(f)(z)  -  a)(z)  -  (z  -  z)$(z) 


F  (am  ,am  )  =  Img 
2  XX  yy  ^ 


n(J)(z)  -  u(z)  -  (z  -  z)$(z) 


The  displacements,  Fj  and  F^,  in  the  composite  sheet  due 
to  the  remote  stress  were  found  as  follows.  First,  the  constitutive 
equations  for  an  orthotropic  material  were  used  to  relate  strainsto 
stresses  (Lekhnitskii  1968)  as 


oc  V 

XX  yx 


c  ^y  E  ^  e 

V  X  y 


where  E  and  E  are  the  moduli  of  elasticity  in  the  x  and  y 
X  y 

directions  respectively,  and  v  and  v  are  Poisson's  ratios. 

./A 

Then  with  the  definition  of  strains  as  ~  lx  ^y  "  W 

equations  (15)  were  integrated  to  find  displacements  as 


/^^xx  ^yx  ( 

u  =  ~  ™yyj  * 


(  X  y 


xy  yy  / 

—  + - >y  +  h2(x) 


32 


where  hl(y)  and  h2(x)  are  arbitrary  functions  which  were  set  to 
zero  because  of  symmetry  considerations. 

The  displacements,  and  F^,  in  the  metal  adherend 
caused  by  the  interlaminar  stresses  were  calculated  by  assuming 
that  the  interlaminar  stresses  acted  as  body  forces  on  the  adherends. 
With  this  assumption,  the  displacements  were  calculated  using  Green's 
functions  in  surface  integrals  as 


Fi(fxz’^yz5  =  aa^x  ""xy 

^Z^^XZ’^yz^  -  ®®yy 


(17) 


where 


“xx  =  "s  '5‘>xx(^>^o>'-fx2<^o”''’'o‘‘^o 

(18) 

=  r/j  SDjy(a.z„)[-fj,^{Zo)]<)x„dy„ 

and  Zg  is  the  location  of  a  point  load  (see  Appendix  D)  and  GD^^, 

GD  .  GD  ,  and  GD  are  the  Green's  functions  which  were 
xy’  yx’  yy 

discussed  and  derived  in  Appendix  D  and  found  as 


GDxx  =  CQ{Real(g(z,ZQ)  +  glz/z^))} 

GDxy  =  CQ{Real(i{g(z,ZQ)  -  g(z,ZQ)})} 

GDyx  =  CQ{Img{g(z,ZQ)  +  g(z,ZQ))} 

"  CQ{Img(i(g(z,ZQ)  -  gCz.z^j))} 

glz.Zp)  =  n[Real{XA(z,ZQ)}  -  XC(z,Zq)] 

+  0.5[n^XB(zjQ)  +  XBIzJq)]  +  XCU.Zq) 

^^0  ■  V 

+  2  — -  +  (z  _  I)B(z,z  ) 

z  - 
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F  (f  ,f  )  =  bb  +  bb 
7^  xz’  yz^  XX  xy 


F  (f  ,f  )  =  bb  +  bb 
xz’  yz^  yx  yy 


where 


‘">xx  =  «'>xx<^k-''k>fxz<"k>'*Vyc 


“vx  =  ™yx<^k""k>fxz<"k>'’'‘o<*>'( 


bbyj,  =  //  HDyy(Z|^,W|^)f^2(W|^)dx„<ly„ 


k  =  1,2 


and  HD  ,  HD  ,  HD  ,  and  HD  are  the  Green's  functionsfor  the 
X  X  xy  y  ^  yy 

composite  adherend  derived  in  Appendix  D  and  found  to  be 


HD^^  =  2Rea1{p  Cllg  (z  ,wj.+  P2C21gj (z^.Wa)} 

AA  ^  1  '*■ 


HD^y  =  2Real{i[PiCl2g2{Zj,wJ  +  P2C22g2(z,  ,w,)]} 


HD^^  =  2Real{qiCngi(z,.Wj)  +  q2C21gi{z^_,v;2)} 

yA 


HDyy  =  2Real[i[qiC21g2(Zi>Wi)  +  q2C22g2(z2,w  J]] 
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z,  =  X  +  s^y,  =  X  +  s,y,  +  s^y^, 


where  s^  and  s^  are  complex  numbers  which  are  not  complex  conjugates 
of  each  other  and  are  roots  of  the  equation  (Lekhnitski  1968) 


Ex  (  E, 

-  -  *  I  - 

’xy  J  Ey 


Cll  = 


s  (s,^v  +1) 

1'  2  '^xy 

4TTtc(s^2  -  $2=^) 


S  +  V 

2  xy 


4TTt  (Sj2  -  ^2  f 


s  is  +  1 ) 

1^  2  VX  ^ 


C21  = 


4Trtc(s^2  - 


S  +  V. 


C22  =  - 


4Trtc{Si^  -  S2^) 


P.  '  r 


Px  =  -  (  =  -  V' 

^x 


1  y  ^  y 
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Substituting  equations  (14),  (16),  (17),  and  (20)  into  equation  (9) 
and  that  result  into  equation  (8)  yields  the  governing  equation  that 
was  used  to  formulate  a  system  of  simultaneous  equations. 

Numerical  Solution 

To  solve  for  the  unknown  stresses  f  and  f  ,  the  domain  of 

XZ  Jr  Z 

integration  in  equations  (18)  and  (21)  was  separated  into  three 
regions  as  shown  on  figure  8.  Region  A  on  the  figure  represents 
the  portion  of  the  domain  where  debonding  has  occurred.  In  this  region 
the  interlaminar  stresses  are  zero.  Region  B  on  the  figure  represents 
the  portion  of  the  domain  where  significant  interlaminar  stresses 
exist.  As  shown  on  the  figure,  this  region  is  divided  into  smaller 
elements.  Region  C  on  the  figure  represents  a  portion  of  the  domain 
where  the  interlaminar  stresses  are  small  and  can  be  neglected.  The 
size  of  each  of  these  regions  will  be  discussed  further  in  Chapter  VI 
where  convergence  of  the  system  is  investigated. 

The  next  step  in  the  formulation  of  the  simultaneous  equations 
was  to  assume  that  the  interlaminar  stresses  were  constant  over  each 
element  of  region  B.  With  this  assumption,  the  displacements  caused 
by  the  interlaminar  stresses,  equations  (17)  and  (20),  were  written 
as 
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y 


Fig.  8.  Domain  of  integration 
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n 

F^(i)  =  Z  AA^^(iJ)[-f^^(j)]  +  AA^y(i,j)[-fy2(j)] 
j=l 


Fg(i)  =  Z  AA^^(iJ)[-fy2(j)]  +  AAyy(i,j)[-fy2(j)] 

j=l 


(23) 


n 

F,(i)  =  Z  BB^y(iJ)[-f^2(j)]  +  BB^y(i,j)[fy2(j)] 

j=l 


n 

F^(i)  =  Z  BBy^(i,j)[f^2(j)]  +  BByy(i,j)[fy2(j)] 
j=1 


where  n  is  the  number  of  elements  in  region  B  and 


AA^^(i,J)  =  /;  GD^^(z,,z„)dx„dy„ 
AA^yd.j)  =  !S  GD^y(z,,z„)dXj,dy^ 

AAy^dJ)  =  //  GD^^(z,,Zj,)dx^dy„ 
AA^yd.j)  =  IS  GDyj,(z,,z„)dXj,dy^ 
BB^^(i,j)  =  II  HD^^(z^,.'\)dx„dy„ 
BB^yd.j)  =  IS  ™^y(2(,i.«|,)dx„dy^ 
BBy^d.j)  =  II  HDj,^(Z|^,,W|^)dx^dy„ 


BB^^(i,j)  =  //  HDy^(Z|.,  ,W|^)dx„dy„ 


k  =  1,2 
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where  the  i  subscript  indicates  the  point  in  the  x-y  plane,  z, 

where  the  governing  equations  (equations  (8))  were  evaluated.  The 

j  subscript  when  used  in  the  coefficient  of  the  interlaminar  stresses, 

f  and  f  ,  indicates  the  location  of  the  centroid  of  the  element 
xz  yz’ 

over  which  the  interlaminar  stresses  act  while  the  j  subscript 
used  in  the  interlaminar  stresses  indicates  the  value  of  these  stresses 
acting  on  element  j.  The  integrals  in  equations  (24)  were  evaluated 
numerically;  the  method  of  integration  will  be  discussed  in  Chapter  VI. 

Substituting  equation  (14)  and  (23)  in  equations  (9)  and  (8)  and 
evaluating  the  result  at  the  centroid  of  each  element  of  region 
B  lead  to  a  system  of  2n  simultaneous  linear  equations  where  n 
is  the  number  of  elements  in  region  B  as 


AA^x(^’j)  +  BB^x(i,j) 

AAxy(lJ)  +  BB^y(' 

i.j) 

+  BBy^( 

ij) 

“yy(’-j)  *  ®®yy' 

ij) 

^ad 

1  0 

F  (i)  -  F  (i) 

1  3 

+  - 

r: 

^eff 

0  1 

_F^(i)  -  F,(i) 

Using  Gaussian  elimination,  this  system  of  linear  simultaneous 
equations  was  solved  and  yielded  values  of  the  unknown  interlaminar 

stresses  ^yz^^^* 
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Nonlinear  Solution 

As  shown  in  Appendix  B,  the  adhesive  used  in  the  reinforced 
system  can  exhibit  nonlinear  stress  strain  behavior.  As  shown  in  the 
appendix,  the  tensile  stress  strain  curve  for  the  adhesive  can  be 
approximated  by  a  bilinear  stress  strain  curve  with  a  change  of  slope 
occurring  at  about  4200  psi.  Because  the  adhesive  in  the  reinforced 
system  is  assumed  to  be  in  a  state  of  pure  shear,  the  data  from  the 
tensile  stress  strain  curve  must  be  related  to  the  adhesive  in  a 
pure  shear  state.  To  develop  this  relationship,  a  yielding  criterion 
is  required. 

For  simplicity  the  criterion  developed  by  Von  Mises  and  given 
by  Hill  (1951)  as 


<“xx  -  °yy>’  <“yy  '  ^  '  "xx'^ 


(26) 


where  a  is  the  stress  at  yielding  and  is  a  constant,  was  used  to 
estimate  when  the  adhesive  in  the  reinforced  system  would  exhibit 
nonlinear  behavior.^  For  pure  tension,  as  was  the  case  in  the  bulk 
property  test  described  in  Appendix  B,  equation  (26)  reduces  to 


^Several  yield  criteria,  of  which  Von  Mises'  is  one  of  the  most 
popular,  are  available  in  the  literature. 
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(27) 


For  the  adhesive  in  the  reinforced  system  the  only  stresses  present, 
according  to  the  Volkerson  assumptions,  are  and  Hence,  in 

this  case  equation  (26)  reduces  to 


(a^  + 

^  yz 


=  V 


(28) 


Combining  equations  (27)  and  (28)  by  eliminating  yields  a  relation 
between  tensile  and  shear  yielding  as 


a 


2 

yz 


+  a" 


xz 


(29) 


For  the  bulk  property  tensile  tests  was  found  approximately 

at  4200  psi.  With  this  value  in  equation  (29)  and  the  notation  for 
interlaminar  stresses  in  the  reinforced  system  an  inequality  was 
developed  as 


f2  +  f2  >  5.88  X  10^  (30) 

yz  xz  - 

Equation  (30)  was  used  to  estimate  the  initiation  of  nonlinear  behavior 
of  the  adhesive  in  the  reinforced  system  by  using  fyz  ^xz 

from  the  solution  of  equation  (25).  As  will  be  shown  in  Chapter  VI  I, 
equation  (30)  predicts  that  nonlinear  behavior  of  the  adhesive  will 
occur  in  many  reinforced  systems.  The  adhesive  nonlinearity  manifests 
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itself  in  equation  (25)  as  changes  in  6^^^  with  the  magnitude  of 
applied  load.  As  shown  in  Appendix  B,  the  for  the  adhesive 

can  be  either  65,000  psi  or  36,000  psi.  The  value  takes  in  the 

reinforced  system  can  be  determined  from  equation  (30).  If  equation  (30) 
is  true  then  ^0ff  “  36,000  psi  while  if  equation  (30)  is  false 
Gg^^  =  65,000  psi. 

By  the  use  of  inequality  (30)  and  equation  (25),  the  applied 
stress  at  which  G^^^  of  the  adhesive  changes  in  the  reinforced 
system  was  predicted  with  the  following  approach.  As  shown  by  equations 
(14)  and  (16)  the  right  hand  side  of  equation  (25)  is  a  function  of 
the  applied  stress,  s  (c^nixx’  ‘^'"yy’  "^^yy  linear  functions 

of  s).  Hence,  the  solution  of  equation  (25)  was  written  in  terms 
of  the  unit  solution  and  the  applied  stress,  s,  as 

fxz(3)  =  ^^xz^-^^unit’  “  ^'^'yz^'^^unit  (31) 

where  s  is  the  remote  applied  stress  and  f^z^^^unit  ^yz^^^unit 

are  the  solutions  of  equation  (25)  with  an  applied  stress  of  s  =  1. 
Substituting  equations  (31)  into  equation  (30)  and  solving  for  s 
for  each  of  the  elements  of  region  B,  lead  to  values  of  the  remote 
stress  s(j)  which  cause  a  change  in  G^^^  for  element  j  as 

5.88  X  10® 


s(j)  = 


(fxz^j^unit^^'*’  ^^z^'^^unit^ 


(37) 
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The  smallest  value  of  s(j)  for  n  elements  of  region  B  is  the  value 
of  the  applied  stress,  s^^^,  which  causes  a  change  in  in 

element  k.  At  this  value  of  s^^^  all  elements  except  element  k 
have  =  65,000,  while  element  k  has  =  36,000. 

If  the  value  of  s^^^  is  greater  than  the  maximum  remote 
applied  stress,  s  ,  then  the  solution  is  completely  elastic  and 

IIIQ  A 

from  equation  (31)  the  interlaminar  stresses  in  the  system  are  given 
as 

fxz(j)  "  ^max^xz^"^  ^unit 

(33) 

fy2(j)  “  ^max^yz^’^^unit 

However,  if  s^^^  is  less  than  s^^^^  then  yielding  has  occurred.^ 

At  the  yield  point  the  stress  in  all  of  the  elements  is  given  by 

9x2  ^  ^xz^'^^unit(i) 

(34) 

where  the  unit(i)  indicates  that  the  unit  stresses  were  obtained  from 
an  elastic  solution  where  all  had  the  same  shear  modulus  of  0^^^ 

Yielding  refers  to  a  change  in  G^^^  in  the  adhesive  of  the 
reinforced  system. 
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65,000  psi.  Once  yielding  has  occurred  in  element  k  the  shear 
modulus  of  element  k  takes  on  the  secondary  value  6^^^  =  36,000. 

Consequently,  for  the  next  increment  of  applied  stress  the 
governing  equations  (25)  must  be  modified  as 


where 


[A]  = 

AAy^(i,j)  +  BByj^(i,j) 


AA^y(i,j)  +  BB^y(i,j) 
AAyy(i  J)  +  BByy(i,j) 


Equation  (35)  was  then  used  to  find  a  new  unit  solution  after 

element  k  had  yielded.  The  new  unit  values,  fxz^^^unit(2) 

f  (i)  were  then  used  in  equation  (31)  and  added  to  equations 

yz'  'umt(2) 

(34)  to  give  the  stress  in  each  element  after  the  second  load  increment 
as 


f„(j) 


^xz^'^^unit(2) 


(36) 
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Substituting  equations  (36)  into  equation  (30)  and  solving  for  s 
yields  the  value  of  s^^^  that  causesthe  next  element  to  yield 


where 


(2)  "^0  V  ^0  "  ^^0^0 


^0  ~  V  xz^'^^unit(2/  {^yz^'^^unit(2) 


-  0^  /3 

yy 


If  the  value  of  s^^^  >  s  then  the  stresses  after  the  section 


increment  of  load  are  given  by 


(Om 


9xz  >  fxz(J)un1t(2)  *  9xz  '<9' 


^vz  ~  ^^max  '  ^yz^'^^unit(2)  '''  ^yz 
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(  2  ) 

and  the  total  number  of  load  increments  is  2.  But  if  s^  '  is  less 

than  s  then  the  stresses  after  the  second  load  increment  are 
max 

given  by 


The  entire  process  is  repeated  until  the  sum  of  all  the  load  increments 

equals  s  .  The  final  value  of  the  interlaminar  stresses  are  the  final 
^  max 

values  of  and 

f,z(j)  =  gxz''^(j) 

(40) 

fyz(j)  =  gyz"’^(j) 

where  m  is  the  total  number  of  load  increments. 

To  obtain  the  final  values  for  the  interlaminar  stresses,  equation 

(40),  the  system  of  simultaneous  equations  (equations  (35))  must  be 

solved  for  each  load  increment.  To  minimize  computational  effort,  a 

Gauss  Siedel  method  discussed  by  McCracken  and  Dorn  (1968)  was  used  to 

solve  equation  (35)  after  the  first  unit  vector  was  found  by  Gaussian 

elimination.  By  the  use  of  unit  vectors,  ^xz^'^^unit(k) 

f  (j)  of  the  k  iteration  as  initial  estimates  for  the  k  +  1 

yz  'um  t(  k ) 

unit  vectors,  the  Gauss  Siedel  method  rapidly  converges.  Consequently, 

the  method  is  very  efficient  in  solving  for  the  unit  stresses  in 
successive  load  increments. 
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In  this  chapter  an  analysis  was  developed  to  determine  inter¬ 
laminar  stresses  in  the  reinforced  system.  In  the  next  chapter  these 
interlaminar  stresses  will  be  used  to  determine  crack  and  debond 
growth  rates  in  the  reinforced  system.  The  accuracy  of  the  analysis 
will  be  discussed  in  Chapter  VI. 


CHAPTER  V 


FATIGUE  ANALYSIS 

As  shown  in  Chapter  II,  fatigue  damage  in  the  reinforced  system 
occurred  as  crack  and  debond  growth.  Accordingly,  an  adequate  fatigue 
analysis  should  predict  both  types  of  growth.  To  the  author's 
knowledge,  the  only  analysis  to  date  that  attempts  to  model  this  behavior 
was  developed  by  Ratwani  (1977).  Ratwani  (1977)  analyzes  crack  and 
debond  growth  in  a  cracked  metal  sheet  reinforced  with  an  uncracked 
metal  sheet  by  using  the  elastic  analysis  developed  by  Erdogan  and 
Arin  (1972)  and  a  maximum  strain  criterion  to  predict  debond  growth. 

In  contrast,  herein,  the  nonlinear  elastic  stress  analysis  developed 
in  the  previous  chapter  was  coupled  with  a  debond  propagation  equation 
based  on  calculated  strain  energy  release  rates. 

Crack  Growth  Rate 

To  use  the  crack  propagation  equation  (equation  (4)),  the  stress 
intensity  must  be  determined.  The  stress  intensity  can  be  related  to 
the  two  stress  distributions  discussed  in  Chapter  IV:  the  remote 
inplane  stresses  in  the  reinforced  system  and  the  interlaminar  shear 
stresses  which  act  near  the  debond  front.  The  stress  intensity  can  be 
found  by  superimposing  the  stress  intensities  from  the  two  stress 
distributions. 
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The  stress  intensity  produced  by  the  remote  stresses  can  be 
calculated  by  substituting  the  stress  function,  $(z),  for  this 
loading  case  (equation  (13))  into  equation  (1)  which  related  the 
stress  function  to  the  stress  intensity  as 


which  results  in 


(42) 


k  =  0 
2 

The  stress  intensity  produced  by  the  interlaminar  shearing 
stresses,  f^_(j)  and  was  found  in  the  following  manner. 

First,  the  stress  intensity  for  four  point  loads  acting  on  a  cracked 
sheet  was  found  by  substituting  the  stress  function 


$(z,Zq)  =  SB(z,Zq)  +  SB(z,Zq) 


derived  in  Appendix  D  into  equation  (1)  to  yield 


(D.43) 
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ki  -  ikj  =  2/2 


lim 

z  ->  a 


/z 


SB(z,Zq)  +  SB(z,Zq) 


(43) 


Taking  the  limit  of  equation  (43)  and  combining  coefficients  of  the 
X  and  Y  load  components  results  in 


k^  -  ik^  =  X[XK(Zq)  +  XK(Zq)] 

+  iY[XK(zQ)  -  XKd^)] 


which  leads  to 


kj  =  2Real[XK(zQ)]X  -  2Img[XK(zQ)]Y 

kj  =  0 


(44) 


where 


/a 


XK(zJ  = 


it(1  +n)t^ 


rz„z„  -  zzJ^  +  a 
0  0  0 

,  (Zo“  - 


nl(Zo)  1 


(45) 


Then,  with  the  use  of  the  coefficient  of  the  X  and 
components  as  Green's  functions  for  the  stress  intensities, 


Y  load 
the  stress 
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intensity  caused  by  the  interlaminar  stress, 
can  be  found  as 


n 

kj  =  2  E  XK(Zq)  dxjjdy^ 

i=l 


+  fy2(j)-^/I"’g(XK(zQ)dXQdy^] 


(46) 


The  domain  of  integration  for  equation  (46)  is,  of  course,  the  same  as 
was  discussed  in  the  previous  chapter.  The  method  of  integration  will 
be  discussed  in  the  next  chapter. 

Debond  Growth  Rate 

To  use  the  debond  propagation  equation  (equation  (5)),  the 
strain  energy  release  rate  must  be  determined  for  points  along  the 
debond  front.  A  rigorous  determination  of  the  strain  energy  release 
rate  is  difficult  and  beyond  the  scope  of  this  effort.  However,  an 
approximate  solution  is  developed  in  the  following  paragraphs. 

A  freebody  of  a  strip  was  taken  from  the  longitudinal  centerline 
of  the  reinforced  system  as  shown  in  figure  9.  The  strain  energy 
release  rate  for  the  strip  was  approximately  calculated  with  the  J 
integral  developed  by  Rice  and  given  in  Liebowitz  (1968)  as 

a  =  J  = 


(47) 


Freebody  for  determination  of  G  along  debond  front 
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where  a.,  is  the  stress,  e . •  is  the  strain,  is  the  surface 
traction,  u^.  is  the  displacement  and  r  is  any  contour  that 
contains  the  debond  front  and  does  not  pass  through  a  plastic  region 
of  material.  Equation  (47)  can  be  written  in  terms  of  cartesian 
coordinates  as  (Yoder  and  Griffis  1974) 


where  is  the  strain  energy  density  and  v  and  w  are  displacements 
in  y  and  z  directions  respectively. 

To  apply  integral  (47)  to  the  freebody  shown  in  figure  9,  a  path 
of  integration  shown  in  figure  10  was  used  to  analyze  the  energy 
release  rate.  The  path  surrounds  the  debond  front  on  which  no 
stresses  or  traction  act,  follows  the  bond  line  in  the  metal  adherend 
on  which  the  interlaminar  stresses  act,  crosses  the  adhesive  which  is 
assumed  to  have  insigificant  stresses,  and  follows  the  bond  line  in  the 
composite  adherend  on  which  the  interlaminar  stresses  act.  The 
SV„ 

-  is  the  strain  in  the  metal  adherend,  and  the  —  is  the 

ay  9y 


strain  in  the  composite  adherend.  With  the  use  of  this  contour,  the 
value  of  G  from  equation  (48)  can  be  written  in  terms  of  the  inter¬ 
laminar  stresses,  f  _  and  f  ,  and  the  strain  in  the  adherends, 


as 


crack 


Fig.  10.  Contour  for  determination  of  G  for  strip 
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With  the  use  of  discrete  values  of  the  interlaminar  stresses  and  strains 
in  the  adherends,  the  strain  energy  release  rate  G  can  be  approximated 
as 

P 

G=  I  Vi)fy,(i)^/i)  -  c„^(i))  (50) 

i=l 

where  represents  the  length  of  the  discrete  elements,  i  the  index 
for  the  elements  in  the  strip,  and  p  the  number  of  elements  in  the 
strip. 

The  interlaminar  stresses,  were  determined  from 

equation  (40).  However,  the  strains  in  the  adherends  still  need  to 
be  determined.  The  strains  in  the  adherends  were  determined  from  the 
stresses  in  the  adherends  according  to  laminate  theory  as  shown  in 
Appendix  C  as 


metal  composite 


(51) 


57 


where  matrices  C  and  D  are  given  by  equations  (C.4)  and  (C.5) 
of  Appendix  C. 

Before  equation  (51)  was  used  to  calculate  the  strain  in  the 
adherends,  the  stresses  in  the  adherends  were  determined. 

As  shown  by  Muskelishvili  (1975)  the  stresses  in  the  metal 
adherend  can  be  expressed  in  terms  of  two  functions,  ^(z)  and  f2(z), 
as 


0  =  Real{3$(z)  -  U{z)  -  (z  -  z)$' (z)}  (52) 

A 

Oy  =  Real{il>(z)  +  ^2(7)  +  (z  -  z)$'  (z)}  (53) 

a  =  Imag{$(z)  +  ^{z)  +  (z  -  z)4>'(z)}  (54) 

xy 

Equations  (52)  -  (54)  were  used  to  determine  the  stresses  in  the  metal 
adherends  caused  by  both  inplane  remote  stresses  and  interlaminar 
stresses. 

For  remote  inplane  stresses,  $(z)  is  given  by  equation  (13) 
as 


4>'(z)  =  - 


(55) 
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and  a{z)  was  found  by  differentiating  equation  (12)  as 


For  the  interlaminar  stresses,  $(z)  and  '}'{z)  were  determined 
and  substituted  in  equations  (52)  -  (54)  (see  Appendix  D).  The  result, 
stress  in  the  metal  caused  by  interlaminar  stresses,  was  found  as 


j=l 


n 


j  =  l 


(57) 

(58) 

(59) 


where  GS^^,  GS^^,  GS^^,  GS^^,  and  are  the  Green's 

functions  for  stresses  given  in  Appendix  D  by  equations  (D.68) 
through  (D.73)  respectively.  The  domain  of  integration  in  equations 
(57)  -  (59)  is  the  area  of  an  element  shown  in  region  B  of  figure  8 


where  j  denotes  the  particular  element. 

Adding  equations  (52)  -  (52)  and  (57)  -  (59)  yielded  the  stress 
in  the  metal  adherend  of  the  reinforced  system.  This  result  was  then 


used  in  equation  (51)  to  calculate  the  strain  in  the  metal  adherend. 


r 
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The  inplane  stresses  in  the  orthotropic  adherend  caused  by  the 
remote  stresses  were  calculated  with  simple  laminate  theory  as  shown 
in  Appendix  C. 

The  inplane  stresses  in  the  orthotropic  adherend  caused  by  the 
interlaminar  stresses  were  found  in  Appendix  F  as 

n 

n 

j=l 
n 

®xy  ”  ^  •^•^^^^(xy)x^xz^'^^  *  *^^(xy)y^yz^'^^^^^o‘^^o 

j=l 

where  HS^^.  HS^^,  HS^^.  and  are  the  Green's 

functions  for  an  orthotropic  solid  sheet  given  in  Appendix  E  by 
equations  (E.23)  -  (E.28). 

Adding  the  inplane  stresses  in  the  orthotropic  adherend  caused 
by  the  remote  stresses  and  the  interlaminar  stresses  and  substituting 
the  result  into  equation  (51)  yielded  the  strains  in  the  orthotropic 
adherend. 

With  the  strains  in  the  adherends,  equation  (50)  was  used  to 
calculate  the  strain  energy  release  rate  for  the  debond  along  the 
longitudinal  axis  of  the  reinforced  system. 


(60) 


(61) 


(62) 
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As  indicated  in  Chapter  III,  the  shape  of  the  debond  throughout 
the  cyclic  tests  can  be  approximated  by  an  ellipse.  The  equation 
for  a  general  ellipse  is  given  as 


y  =  b[i  -  (f)^] 


ih) 


(63) 


where  a  is  half  the  crack  length  and  b  is  the  debond  length 
measured  along  the  y-axis  from  the  center  of  the  crack.  With  the  use 
of  equation  (63)  the  overall  debond  shape  was  predicted  by  calculating 
the  half  crack  length,  a,  from  equation  (4)  using  the  stress 
intensity  calculated  from  equation  (45)  and  calculating  the  debond 
length,  b,  from  equation  (5)  using  the  strain  energy  release  rate 
calculated  from  equation  (50). 

Prediction  of  Crack  and  Debond  Growth 

The  analysis  developed  in  this  dissertation  was  programmed  on 
a  CDC  6600  computer  at  NASA-LRC.  A  discussion  of  the  program,  a 
sample  analysis,  and  a  program  listing  are  given  in  Appendix  F. 

A  flow  chart  of  key  elements  of  the  program  is  shown  in  figure  11. 

In  the  next  chapter  the  accuracy  of  the  analysis  is  investigated 
by  comparing  results  of  the  calculations  from  the  analysis  with 
experimental  results  generated  in  Chapter  III. 
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Fig.  11.  Flow  chart  for  fatigue  analysis 


CHAPTER  VI 


ASSESSMENT  OF  ACCURACY 

Before  the  analysis  developed  in  the  previous  chapters  was  used 
to  study  the  fatigue  behavior  of  the  reinforced  systems,  the  accuracy 
of  the  analysis  was  assessed.  This  assessment  was  based  on  analytical 
convergence  studies  and  studies  comparing  analytical  results  with 
experimental  results  obtained  in  Chapter  III. 

Numerical  Integration 

A  key  item  in  the  analysis  was  the  method  of  integration  of  the 
Green's  functions,  In  theory  the  Green's  functions  for  both  stress 
and  displacement  require  integration  over  an  infinite  domain.  Because 
the  functions  are  complicated,  a  closed  form  integration  is  difficult 
if  not  impossible  to  perform.  Consequently,  a  numerical  solution  was 
employed.  Two  key  items  in  this  numerical  integration  were  the  domain 
of  integration  and  the  method  of  numerical  integration. 

As  shown  on  figure  8,  the  infinite  domain  of  integration  was 
divided  into  three  regions:  A,  B,  and  C,  Only  region  B  has  significant 
interlaminar  shear  stresses  which  need  be  integrated.  To  perform  the 
integration  region  B  was  divided  into  elements  as  shown  on  figure  8. 

Each  of  these  elements  was  bounded  by  the  curves 
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X  =  Xj  X  =  Xj 

yi  =  fi(x)  ya  =  f2(x) 


(64) 


where  the  functions  fj  and  fz  are  the  form  of  equation  (63)  and 
reflect  the  debond  shape.  For  example,  the  bounding  curves  for 
elements  along  the  debond  front  are  given  as 


yjx)  =  b  /I  - 


y2(x)  =  (b  +  db) 


i’l- 


da 


(65) 


where  a  and  b  are  the  half  crack  length  and  the  debond  height 
respectively  and  da  and  db  represent  fractions  of  a  and  b. 

The  numerical  method  of  integration  used  for  each  element  of 
the  region  B  was  a  two-dimensional  Simpson's  integration.  Simpson's 
integration  was  used  so  that  several  of  the  Green's  functions  given  by 
equations  (19),  (22),  (D.68)  -  (D.73),  and  (E.23)  -  (E.28)  could  be 
integrated  simultaniously  by  using  common  values  of  the  complex 
functions.  In  each  element  the  interlaminar  stresses  were  assumed 
to  be  constant  and  the  Green's  functions  were  integrated  by  using  9  or 
18  integration  points.  Nine  integration  points  were  used  when  the 
domain  of  integration  did  not  contain  a  singularity  while  18  points 
were  used  when  a  singularity  existed  within  the  element. 
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In  the  latter  case  the  domain  of  the  element  was  divided  into 
three  regions;  two  of  which  were  analytic  and  the  third,  which  contained 
the  singularity,  was  nonanalytic.  The  two  analytic  regions  were 
integrated  with  the  nine  point  Simpson's  integration  scheme.  The 
nonanalytic  region  was  integrated  by  separating  the  Green's  functions 
into  products  of  analytic  and  nonanalytic  functions.  The  analytical 
portions  were  expanded  in  a  Taylor  series  and  only  the  first  terms, 
which  were  constants,  were  retained.  The  singular  portion  was  integrated 
analytically  in  the  principal  value  sense  (Hilderbrand  1950).  The 
product  of  the  first  term  of  the  Taylor  series  and  the  principal 
value  resulted  in  an  approximate  integral  value  in  the  nonanalytic 
region.  In  general  the  value  of  the  integral  in  the  nonanalytic  region 
of  the  element  was  small  in  comparison  with  values  in  the  two  analytic 
regions. 

The  size  of  region  B  shown  in  figure  8  was  determined  iteratively 
by  starting  with  a  small  region  and  increasing  its  size  until  no 
changes  occurred  in  the  interlaminar  stresses,  stress  intensity  or 
strain  energy  release  rate.  As  an  example.  Panel  B  (under  a  load  of 
22,500  pounds)  discussed  in  Chapter  III  was  modeled  as  shown  in  figure 
12  with  an  initial  crack  length  of  one  inch  and  no  debond  (b  =  0).  As 
shown  in  the  figure,  region  B  was  assumed  to  be  nearly  rectangular 
with  n  elements  along  the  x-axis  and  n  elements  along  the  y-axis. 

X 

The  elements  have  a  width  of  0.25  inch  and  a  height  of  0.20  inch.  The 
analysis  was  performed  for  values  of  n^  ranging  from  four  to  six  and 
values  of  n^  ranging  from  one  to  six.  Test  conditions  for  each  of 
these  possible  combinations  were  notated  as 
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Fig.  12.  Domain  B  for  convergence  analysis 
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Figure  13  shows  the  values  of  interlaminar  stresses  along  the 

longitudinal  centerline  for  the  different  conditions.  As  evident 

from  the  figure, values  of  the  interlaminar  stresses  have  converged 

for  values  of  n  =  4  and  n  =  3.  Also  both  the  stress  intensities 

X  y 

and  the  strain  energy  release  rates  have  converged  for  these  values. 
Thus,  for  this  sample  case  the  domain  of  region  B  is  about  1  x  0.6  inch. 
Similar  analyses  showed  that  the  length  of  the  domain  of  region  B 
along  the  x-axis  is  typically  the  length  of  the  crack.  However,  as  the 
adherend  thickness  or  shear  modulus  of  the  adhesive  changes, the 
extent  of  region  B  in  the  y-direction  also  changes. 

With  the  use  of  the  one-dimensional  analysis  discussed  in 
Appendix  B,  the  extent  of  region  B  in  the  y-direction  was  estimated  for 
different  adherend  thicknesses  and  adhesive  moduli  in  the  following 
manner.  A  strip  was  taken  from  along  the  longitudinal  centerline  of  the 


interlaminar  shear  stresses,  psi 
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Fig.  13.  Convergence  of  interlaminar  shear  stresses 
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reinforced  system  as  shown  in  figure  9.  Because  of  the  crack  the 
load  in  this  strip  must  be  transferred  either  to  the  adjacent  metal 
via  inplane  shear  stresses  or  to  the  composite  via  interlaminar  stresses. 
If  all  the  load  were  transferred  via  the  interlaminar  stresses,  the 
interlaminar  stresses  could  be  higher  and  extend  over  a  greater  area 
than  if  the  inplane  stresses  were  considered.  Consequently,  the  boundary 
of  the  interlaminar  stress  region  calculated  from  the  one-dimensional 
analysis,  which  only  considers  interlaminar  stresses,  would  be  an 
upper  bound. 

From  the  one-dimensional  analysis  the  shearing  stresses  are 
found  by  equation  (B-12)  as 


T(y)  =  (B.12) 

where 

1  ]  s-ad 

+ - )  K  =  - 

Examination  of  equation  (B.12)  revealed  that  the  T(y)  is  a  maximum  at 
y  =  0.  The  shear  stresses  were  assumed  to  be  negligible  when  they 
were  smaller  than  5  percent  of  the  maximum  value.  In  equation  form, 
the  relationship  between  y  and  95  percent  of  the  maximum  inter¬ 
laminar  shear  stress  was  expressed  as 


Gad  (  1 

t  .  )  t  E 

ad  (.mm 
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the  distance  at  which  the  shear  stresses  are  95  percent  of  their 
maximum  value.  As  evident  from  equation  (66)  the  y  distance  is  a 
function  of  the  thicknesses  and  material  properties  of  the  adherends 
and  the  adhesive.  For  the  adhesive  used  in  the  reinforced  system  of 
Panels  A  and  B  two  effective  shear  moduli  were  found  as  65,000  psi 
prior  to  yielding  and  26,100  psi  after  yielding  (see  Appendix  B). 

With  the  use  of  the  latter  of  these  two  values  to  give  a  conservative 
bounds,  the  y  distance  for  region  B  of  Panel  B  (see  figure  2)  as 
calculated  by  equation  (66)  was  found  as  y  =  0.60  in.  The  estimated 
value  agrees  well  with  convergence  study  results  shown  on  figure  13. 

The  same  logic  applies  to  both  bonded  (b  =  0)  and  debonded  systems 
(b  >  0).  Consequently,  the  domain  of  region  B  used  in  the  integration 
of  the  Green's  functions  was  determined  by  the  crack  length  and 
equation  (66). 

Once  the  domain  of  region  B  was  determined,  the  effect  of  mesh 
refinement  within  region  B  was  investigated.  Figure  14  shows  a 
comparison  of  interlaminar  stresses  calculated  (again  for  Panel  B) 
using  two  different  mesh  sizes  in  region  B.  As  evident  from  the  figure. 
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The  effect  of  mesh  size  on  interlaminar  stresses 


71 


the  interlaminar  stresses  calculated  with  the  fine  mesh  can  deviate 
as  much  as  16  percent  from  stresses  calculated  with  the  coarse  mesh.  The 
stress  intensities  and  strain  energy  release  rates  calculated  using 
the  two  different  meshes  were  found  as 

k  G 

1 

coarse  mesh  0.281  1.29 

fine  mesh  0.280  1.36 

The  values  of  and  G  are  not  nearly  as  sensitive  to  changes  in  mesh 
size  as  the  interlaminar  stresses  were.  The  reason  for  the  insensitivity 
is  most  likely  because  both  k^  and  G  are  obtained  by  integration 
schemes  that  smooth  out  the  effect  of  local  approximations  in  the 
interlaminar  stresses.  Consequently,  the  grid  can  be  rather  coarse 
and  still  accurately  predict  both  k^  and  G.  In  contrast  the  values 
of  the  interlaminar  stresses  require  a  finer  mesh  for  accurate  values. 
Because  k^  and  G  are  used  to  predict  the  fatigue  behavior  of  the 
reinforced  system,  a  relatively  coarse  mesh  was  used  in  the  analysis 
without  loss  of  accuracy. 

Accuracy  of  the  Analysis 

To  ascertain  the  accuracy  of  the  analysis,  calculated  values  of 
stresses  in  the  adherends,  the  stress  intensities,  and  the  crack 
propagation  rates  and  debond  sizes  were  compared  to  experimental 
results  on  Panels  A  and  B  shown  in  Chapter  III.  To  compare  calculated 
and  experimental  values  of  stresses  in  the  adherends  and  stress 


intensities.  Panels  A  and  B  were  modeled  when  the  half  crack  length,  a, 
reached  2  inches  as  shown  on  figure  15.  The  debond  sizes  were  obtained 
from  the  C-scans  of  the  specimens  shown  on  figure  5.  The  height  of  the 
grid  from  the  edge  of  the  debond  was  determined  from  equation  (66). 

The  crosshatched  elements  on  the  figure  indicate  elements  in  which  the 
adhesive  has  changed  modulus  (yielded)  during  each  applied  load  cycle. 
The  effect  of  the  modulus  change  on  and  G  will  be  discussed  in 
the  next  chapter. 

% 

The  values  of  the  crack  growth  rate  for  Panels  A  and  B  were 
calculated  for  a  half-crack  length  of  2  inches,  the  debond  sizes 
observed  in  figure  5,  and  the  applied  loads  shown  on  page  18. 

The  calculated  values  and  the  experimental  values,  from  Chapter  II, 
of  the  crack  propagation  rates  are 

da/dN,  in/cycle 

Panel  calculated  experimental 

A  2.85  X  10"^  1.30  X  10"^ 

B  2.26  X  10"^  3.12  x  lO"^ 

The  difference  between  the  calculated  and  experimental  rates  are 
within  the  scatter  of  the  predicted  rates  for  unreinforced  metal 
sheets  (Figge  and  Newman  1967).  Hence,  the  analysis  accurately 
predicts  the  crack  growth  rate  if  the  crack  length  and  debond  are 
known. 

As  a  further  check,  the  analysis  was  used  to  predict  the  crack 
growth  in  Panels  A  and  B.  Figure  16  shows  the  calculated  and 
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Fig.  15.  Mesh  for  analysis  of  Panels  A  and  B 


Panel  A 
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Fig.  16.  Experimental  and  calculated  crack  lengths  as  a  function  of  applied  load  cycles 
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experimental  half-crack  lengths  plotted  against  the  number  of  applied 
load  cycles.  For  both  panels  the  analysis  gave  a  conservative  prediction 
of  the  crack  growth.  The  calculated  and  experimental  half-crack  lengths 
agree  within  a  factor  of  2.  This  deviation  is  within  the  scatter  of 
crack  length  prediction  of  unreinforced  metals.  Hence,  the  analysis 
appears  to  accurately  predict  the  crack  length  as  a  function  of 
applied  load  cycles. 

On  figure  17  the  calculated  and  experimental  crack  propagation 
rates  were  plotted  against  the  half-crack  length.  The  calculated 
crack  propagation  rates  were  within  a  factor  of  2  of  the  experimental 
rates.  The  largest  error  occurred  in  Panel  B  for  a  half-crack  length 
of  2  inches. 

The  calculated  crack  lengths  and  crack  propagation  rates  shown 
on  figures  16  and  17  are  a  function  of  debond  growth.  On  figure  18 
the  debond  aspect  ratio,  b/a,  was  plotted  against  the  half-crack 
length  for  the  two  panels.  The  symbols  on  the  figure  indicate  values 
of  the  debond  aspect  ratio  obtained  experimentally  (see  figure  10). 

As  evident  on  the  figure,  the  calculations  indicated  that  the  debond 
aspect  ratio  increases  rapidly,  especially  for  Panel  A,  before  the 
half-crack  length  reaches  0.2  inch.  Hence,  the  debond  grows  before 
the  crack  does. 

Of  the  two  panels,  Panel  A  exhibits  the  largest  debond  growth,  and 
at  a  half-crack  length  of  2  inches  the  predicted  debond  aspect  ratio  was 
determined  experimentally.  For  Panel  B,  the  analysis  predicted  a  debond 
aspect  ratio  of  0.45  at  a  half-crack  length  of  2  inches.  In  contrast. 


crack  propagation  rate  da/dN,  inches/cycle 
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Fig.  17.  Experimental  and  calculated  crack  propagation  rates 


experimental 


Debond  aspect  ratio  versus  crack  length 
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the  debond  aspect  ratio  obtained  experimentally  was  about  0.14.  The 
discrepancy  between  the  calculated  and  experimental  values  may  be 
linked  to  the  magnitude  of  strain  energy  release  rate  as  the  debond 
propagates  in  Panel  B.  Values  of  G  in  Panel  B  ranged  from 
0.63  in-lbs/in  at  a  half-crack  length  of  1  inch  to  0.53  in-lb/in  at  a 
half-crack  length  of  2  inches  (in  contrast  G  values  for  Panel  A 
ranged  from  6.0  to  2.45  in-lbs/in).  The  values  of  G  for  Panel  B 
were  below  values  of  G  used  to  determine  the  empirical  constants  in 
equation  (5).  A  few  exploratory  tests  of  the  type  performed  in  Appendix 
A  revealed  that  a  threshold  value  of  G  may  exist.  Below  this  threshold 
value  debonding  does  not  occur.  Although  proving  the  existence  of  a 
threshold  was  beyond  the  scope  of  this  dissertation,  if  it  does 
exist  the  analysis  of  Panel  B  would  predict  a  debond  aspect  ratio  much 
closer  to  the  experimental  value. 

The  comparison  between  calculated  and  experimental  values  of 
predicted  crack  length,  crack  propagation  rate  and  debond  aspect 
ratio  shown  on  figures  16  through  18  showed  that  the  analysis  has 
potential  to  predict  crack  growth  in  metals  reinforced  with  composite 
materials.  However,  a  true  assessment  of  the  accuracy  of  the  analysis 
can  be  made  only  after  a  more  extensive  data  base  is  developed. 


CHAPTER  VII 


PARAMETRIC  STUDIES  ON  CRACK  AND  DEBOND  GROWTH 

To  give  insight  about  crack  and  debond  growth  in  reinforced 
systems,  the  analysis  was  performed  for  reinforced  systems  with 
several  different  adherend  thicknesses,  debond  sizes,  and  crack 
lengths.  The  metal  adherend  thicknesses  studied  were  0.05,  0.10,  and 
0.15  inch;  the  composite  adherend  thicknesses  were  0.025,  0.05,  and 
0.075;  the  half-crack  lengths  were  0.5,  1.0,  1.5  inches;  and  the  aspect 
ratios  of  the  debond  areas  were  0.001,  0.5,  and  1.0.  For  the  various 
combinations  of  adherend  thicknesses,  crack  lengths,  and  debond 
aspect  ratios,  the  stress  intensities,  strain  energy  release  rates,  and 
remote  stress  that  caused  nonlinear  behavior  of  the  adhesive  were 
calculated  and  are  shown  in  Table  2. 

The  first  two  columns  of  the  table  give  the  metal  and  adherend 
thicknesses,  t^  and  t^.  The  third  and  fourth  columns  give  the 
half-crack  length,  a,  and  the  debond  aspect  ratio, b/a.  The  fifth 
column  gives  the  stress  intensity  in  the  metal  adherend  normalized  by 
the  remote  stress,  k/s.  The  sixth  column  gives  the  strain  energy 
release  rate  at  the  debond  front  (along  the  longitudinal  centerline) 
normalized  by  the  remote  stress  squared,  G/s^.  The  last  column  gives 
the  remote  stress  that  would  cause  the  adhesive  layer  in  the  reinforced 
system  to  behave  non! i nearly. 
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TABLE  2 

CALCULATED  VALUES  FOR  PARAMETRIC  STUDY 


m 

(in) 

‘c 

(in) 

a 

(in) 

b/a 

k/s 

(in^) 

G/s^ 

(inVlbs) 

^yielid 

psi 

0.05 

50 

0.025 

0.50 

0.001 

0.171 

12.4E-10 

19,600 

0.50 

0.50 

0.240 

8.2E-10 

25,200 

0.50 

1.00 

0.289 

4.4E-10 

25,600 

1.00 

0.001 

0.170 

12.6E-10 

17,600 

1.00 

0.50 

0.286 

9.5E-10 

23,900 

1.00 

1.00 

0.366 

4.9E-10 

23,600 

1.50 

0.001 

0.171 

12.4E-10 

16,500 

1.50 

0.50 

0.322 

9.5E-10 

23,800 

o.o: 

15 

1.50 

1.00 

0.426 

5. IE-10 

23,700 

0.050 

0.50 

0.001 

0.138 

7.3E-10 

25,000 

0.50 

0.50 

0.190 

5.6E-10 

21,700 

0.50 

1.00 

0.230 

3.0E-10 

29,500 

1.00 

0.001 

0.136 

7.2E-10 

22,800 

1.00 

0.50 

0.221 

6.0E-10 

28,100 

1.00 

1.00 

0.286 

3.3E-10 

27,000 

1.50 

0.001 

0.137 

7.2E-10 

21,100 

1.50 

0.50 

0.264 

6.0E-10 

27,600 

0.0! 

50 

1.50 

1.00 

0.329 

3.3E-10 

27,200 

0.075 

0.50 

0.001 

0.125 

5.7E-10 

28,200 

0.0! 

50 

0.075 

0.50 

0.50 

0.170 

4.8E-10 

33,800 
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TABLE  2-Continued 


m 

(in) 

(in) 

a 

(in) 

b/a 

k/s 

(in"'^) 

G/s^ 

(in‘*/lbs) 

^yield 

psi 

0.( 

)50 

0.( 

375 

0.50 

1.00 

0.206 

2.4E-10 

31,900 

1.00 

0.001 

0.123 

5.6E-10 

26,500 

1 

1.00 

0.50 

0.196 

4.7E-10 

30,600 

1.00 

1.00 

0.254 

2.6E-10 

29,100 

1.50 

0.001 

0.123 

5.6E-10 

23,800 

1.50 

0.50 

0.217 

4.8E-10 

30,000 

0.( 

)50 

0.( 

375 

1.50 

1.00 

0.291 

2.7E-10 

29,300 

0.100 

0.( 

325 

0.50 

0.001 

0.254 

41.5E-10 

11,200 

0.50 

0.50 

0.334 

24.6E-10 

15,700 

0.50 

1.00 

0.385 

11.4E-10 

17,600 

1.00 

0.001 

0.256 

46.1E-10 

9,800 

1.00 

0.50 

0.407 

28.9E-10 

14,200 

1.00 

1.00 

0.498 

13.5E-10 

15,400 

1.50 

0.001 

0.257 

46.2E-10 

9,100 

1.50 

0.50 

0.462 

29.9E-10 

13,700 

0.( 

325 

1.50 

1.00 

0.583 

14.2E-10 

15,100 

0.050 

0.50 

0.001 

0.202 

23.9E-10 

14,600 

0.50 

0.50 

0.263 

16.2E-10 

18,900 

0.50 

1.00 

0.307 

8.1E-10 

19,900 

1.00 

0.001 

0.202 

25.1E-10 

13,100 

0.' 

100 

0.( 

350 

1.00 

0.50 

0.309 

18.5E-10 

17,400 
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TABLE  2-Continued 


■'m 

(in) 

0.100 


0.100 

0.150 


(in) 

a 

(in) 

b/a 

k/s 

(in"^^) 

G/s^ 

(in‘‘/lbs) 

^yield 

psi 

0.050 

1.00 

1.00 

0.385 

9.4E-10 

12,379 

1.50 

0.001 

0.201 

25.0E-10 

16,800 

1.50 

0.50 

0.345 

19.0E-10 

16,800 

O.C 

)50 

1.50 

1.00 

0.444 

9.9E-10 

16,700 

0.075 

0.50 

0.001 

0.178 

17.6E-10 

16,900 

0.50 

0.50 

0.230 

12.5E-10 

21,300 

0.50 

1.00 

0.269 

6.4E-10 

21,800 

1.00 

0.001 

0.177 

18.0E-10 

15,200 

1.00 

0.50 

0.266 

14.0E-10 

19,400 

1.00 

1.00 

0.333 

7.4E-10 

18,500 

1.50 

0.001 

0.176 

17.9E-10 

14,600 

1.50 

0.50 

0.294 

14.4E-10 

18,600 

0.( 

175 

1.50 

1.00 

0.382 

7.4E-10 

18,100 

0.025 

0.50 

0.001 

0.313 

78.3E-10 

8,300 

0.50 

0.50 

0.397 

40.9E-10 

12,400 

0.50 

1.00 

0.446 

18.0E-10 

14,800 

1.00 

0.001 

0.323 

95.6E-10 

7,000 

1.00 

0.50 

0.493 

51.3E-10 

10,900 

1.00 

1.00 

0.584 

22.3E-10 

12,500 

1.50 

0.001 

0.325 

99.0E-10 

6,500 

1.50 

0.50 

0.565 

54.1E-10 

10,400 

0.1 

325 

1.50 

1.00 

0.689 

23.7E-10 

12,100 

0.150 
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TABLE  2-Continued 


(in) 

r 

C 

M 

n) 

a 

(in) 

b/a 

k/s 

(in"'^) 

G/s^ 

(in^/lbs) 

^yield 

psi 

0.1 

50 

0.050 

0.50 

0.001 

0.253 

46.8E-10 

10,600 

0.50 

0.50 

0.318 

28.9E-10 

14,400 

0.50 

1.00 

0.362 

13.9E-10 

16,200 

1.00 

0.001 

0.256 

52.9E-10 

9,400 

1.00 

0.50 

0.379 

35.2E-10 

12,900 

1.00 

1.00 

0.460 

17.1E-10 

13,400 

1.50 

0.001 

0.256 

53.7E-10 

8,700 

1.50 

0.50 

0.425 

37.1E-10 

12,400 

0.1 

150 

0.( 

)50 

1.50 

1.00 

0.534 

18.2E-10 

12,900 

0.150 

0./ 

^5 

0.50 

0.001 

0.222 

34.4E-10 

12,300 

0.50 

0.50 

0.278 

22.5E-10 

16,100 

0.50 

1.00 

0.318 

11.3E-10 

17,600 

1.00 

0-001 

0.223 

37.3E-10 

11,100 

1.00 

0.50 

0.325 

26.9E-10 

14,500 

1.00 

1.00 

0.397 

13.7E-10 

14,400 

1.50 

0.001 

0.222 

37.6E-10 

10,300 

1.50 

0.50 

0.360 

28.2E-10 

14,000 

0.' 

150 

0.: 

75 

1.50 

1.00 

0.457 

14.5E-10 

13,700 
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Stress  Intensity  Factor,  Crack  Growth  Rate 

Figures  19,  20,  and  21  show  the  effects  of  composite  adherend 
thickness,  debond  aspect  ratio,  and  crack  lengths  on  stress  intensity 
for  the  different  thickness  metal  adherends.  On  the  figures  circles, 
diamonds  and  squares  represent  stress  intensity  values  for  0.025,  0.05, 
and  0.075  inch  thick  composite  adherends  respectively.  Open,  half- 
shaded,  and  shaded  symbols  represent  stress  intensity  values  for 
debond  aspect  ratios  of  0.0,  0.50,  and  1 .0  respectively.  On  the  figures 
the  stress  intensity  normalized  by  the  remote  stress  was  plotted  against 
the  half-crack  length. 

For  all  three  metal  adherend  thicknesses,  the  figures  show 
consistent  trends.  If  the  debond  aspect  ratio  is  small,  the  stress 
intensity  is  not  significantly  affected  by  either  the  crack  length  or 
the  thickness  of  the  composite  reinforcement.  However,  as  the 
debond  size  increases,  the  stress  intensities  increase  significantly 
with  longer  crack  lengths  and  thinner  composite  reinforcement. 

The  effects  of  debond  size,  composite  adherend  thickness,  and 
crack  length  become  more  pronounced  for  the  thicker  metal  adherends. 

Strain  Energy  Release  Rate,  Debond  Propagation 

Figures  22,  23,  and  24  show  the  effect  of  composite  adherend 
thickness,  debond  aspect  ratio,  and  crack  length  on  the  strain  energy 
release  rate  at  the  debond  front  (along  the  longitudinal  centerline 
of  the  reinforced  system).  As  on  figures  19,  20,  and  21,  different 
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symbols  indicate  different  composite  adherend  thicknesses,  and 
different  amounts  of  shading  indicate  different  aspect  ratios.  On 
the  figures  the  strain  energy  release  rate  normalized  by  the  remote 
stress  squared  was  plotted  against  the  half  crack  length. 

Examination  of  the  figures  revealed  that  for  all  metal  adherend 
thicknesses  the  thickness  of  the  composite  had  the  most  pronounced 
effect  on  the  strain  energy  release  rate.  The  thinner  the  composite 
adherend  the  higher  the  strain  energy  release  rate  and  the  more 
likely  debonding  would  occur.  The  debond  aspect  ratio  had  the  next 
most  significant  effect.  The  larger  the  debond  aspect  ratio  the  lower 
the  strain  energy  release  rate.  Of  all  the  parameters  the  crack 
length  had  the  smallest  effect  on  the  strain  energy  release  rate. 

As  in  the  case  for  stress  intensities,  the  effects  of  debond 
size,  composite  adherend  thickness,  and  crack  length  are  more 
pronounced  for  the  thicker  metal  adherends.  In  fact  on  figure  24 
the  energy  release  rates  for  0.025  inch  composite  reinforcement  with  no 
debond  reinforcing  a  0.15  inch  metal  were  so  great  for  all  crack 
lengths  that  the  energy  release  rates  were  off  scale  in  the  figure. 
Evidently,  the  most  severe  case  for  debond  occurs  with  a  thick 
metal  adherend  reinforced  with  a  thin  composite  sheet  with  no 
debonding  between  adherends. 

Nonlinear  Effects 

Figures  19  through  24  were  generated  by  assuming  that  the  adhesive 
behaved  linearly.  Hence,  because  they  were  normalized  with  respect 
to  the  remote  stress  or  its  square,  they  can  be  used  to  estimate  the 
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stress  intensity  and  strain  energy  release  rate  for  any  remote  stress 
that  does  not  produce  nonlinear  behavior  of  the  adhesive  in  the 
reinforced  system.  For  the  different  parameters  studied,  Table  2 
gives  the  values  for  the  remote  stresses  that  cause  nonlinear  behavior 
of  the  adhesive.  As  evident  from  the  table  nonlinear  behavior  can 
occur  at  relatively  low  remote  stresses.  For  values  in  the  table, 
the  lowest  value  of  remote  stress  to  cause  nonlinear  behavior  occurved 
for  the  thickest  (0.15  in)  metal  adherend  with  a  1.5  inch  crack  that 
was  reinforced  with  the  thinnest  (0.025  in)  composite  adherend  and  no 
debond.  For  this  reinforced  system  a  remote  stress  of  6,000  psi 
caused  nonlinear  behavior  of  the  adhesive.  However,  for  practical 
purposes,  the  remote  stress  applied  to  this  system  could  be  as  high 
as  52,000  psi.  To  investigate  the  effects  of  the  nonlinear  adhesive 
on  the  crack  and  debond  growth  predictions,  the  analysis  was  conducted 
using  the  preceding  parameters  for  both  a  linear  adhesive  and  nonlinear 
adhesive. 

For  the  linear  analysis  an  effective  shear  modulus  (see  Appendix 
B)  of  65,000  psi  was  used,  while  for  the  nonlinear  analysis  an  effective 
shear  modulus  of  65,000  psi  was  used  until  the  remote  si.ress  reached 
6,000  psi  after  which  an  effective  shear  modulus  of  36,000  psi  was 
used  for  nonlinear  elements  of  region  B  (see  page  38).  With  the  use 
of  the  two  different  analyses,  the  stress  intensity  and  the  strain 
energy  release  were  found  as 


93 


adhesive 

G 

linear 

18,800 

26.1 

nonl inear 

16,700 

26.1 

percent  error 

13 

0 

For  this  example,  the  nonlinear  analysis  predicted  the  stress 
intensity  13  percent  below  that  predicted  by  the  linear  analysis. 
However,  the  predicted  strain  energy  release  rate  was  the  same  for 
both  analyses. 

At  first  glance  the  linear  analysis  would  then  accurately 
predict  the  debond  growth  and  conservatively  estimate  the  crack  growth. 
But  by  the  use  of  equations  (4)  and  (5)  and  the  above  values  the  debond 
and  crack  growth  rates  were  found  at  4.2  inches/cycle  and  1.56E-04 
inches/cycle  respectively.  Hence,  the  debond  propagates  much  faster 
than  the  crack.  As  the  debond  grows,  the  magnitude  of  the  stress 
intensity  from  the  linear  and  nonlinear  analysis  converges.  In  fact, 
because  the  debond  grows  so  much  faster  than  the  crack,  before  the 
crack  extends  any  appreciable  amount  the  stress  intensities  from  the 
two  analyses  predict  the  same  crack  growth  rates.  In  addition,  because 
the  example  exhibits  the  most  significant  nonlinearity  of  all  the  cases 
considered  in  Table  2,  the  linear  analysis,  and  hence  figures  19 
through  24  can  be  used  to  estimate  crack  and  debond  growth  even  when 
the  adhesive  behaves  nonlinearly. 
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Prediction  of  Crack  and  Debond  Grov/th 

To  predict  the  crack  and  debond  growth  in  a  reinforced  system, 
the  computer  code  discussed  in  Appendix  F  should  in  general  be  employed. 
However,  figures  19  through  24  can  also  be  used  to  estimate  both  the 
debond  and  crack  growth  for  a  variety  of  reinforced  systems  in  the 
following  manner.  First,  for  the  given  adherend  thickness,  crack 
length,  debond  size,  and  remote  applied  stress,  the  strain  energy- 
release  rate,  G,  and  the  stress  intensity,  k^ ,  can  be  estimated 

from  figures  19  through  24.  Then,  equations  (4)  and  (5),  the  crack 
and  debond  growth  equations,  can  be  used  to  estimate  the  number  of 
applied  load  cycles  to  extend  the  crack  by,  for  exam.ple,  10  percent 
and  to  extend  the  debond  length  (along  the  longitudinal  axis  of  the 
reinforced  system)  by  10  percent.  The  smallest  value  of  the  applied 
load  cycles  required  to  produce  these  extensions  is  used  with  equations 
(4)  and  (5)  to  predict  the  extension  of  crack  and  debond  growth.  Next, 
these  extensions  are  added  to  the  original  crack  and  debond  length 
and  the  entire  process  is  repeated  until  the  stress  intensity  reaches 
56,000  psi-in”^  (fracture  occurs)  or  the  desired  crack  length  or 
number  of  applied  load  cycles  is  reached. 


CONCLUSIONS 


The  failure  mode  of  cracked  metal  sheets  that  are  reinforced 
with  composite  is  crack  propagation  in  the  metal  sheet.  Analysis  of 
the  crack  growth  is  complicated  by  the  development  of  a  debond  near 
the  crack.  Herein,  an  analysis  was  developed  to  predict  both  the  debond 
and  crack  growth  in  a  reinforced  system.  The  analysis  was  predicated 
on  the  use  of  strain  energy  release  rate  to  correlate  debond  growth. 
Empirical  constants  required  for  the  correlation  were  developed 
from  simple  bonded  specimens.  The  correlating  equation  for  the 
debond  growth  was  then  used  in  a  stress  analysis  that  was  based  on 
complex  variable  Green's  functions  which  were  developed  herein  for 
cracked,  isotropic  sheets  and  uncracked,  orthotropic  sheets.  The 
stress  analysis  was  used  to  calculate  the  inplane  and  interlaminar 
stresses,  the  stress  intensity  at  the  crack  tip,  and  the  strain 
energy  release  rate  at  the  debond  front.  By  the  use  of  the  analysis, 
an  iterative  solution  was  developed  that  used  the  stress  intensity  and 
the  strain  energy  release  rate  to  predict  the  crack  and  debond  growth 
on  a  cycle- by- cycle  basis. 

To  verify  the  analysis,  tests  were  conducted  on  two  different 
reinforced  panels  which  exhibited  different  amounts  of  debonding. 

For  both  panels  the  predicted  crack  growth  was  within  the  accuracy 
of  crack  growth  prediction  in  unreinforced  metal  sheets.  Hence, 
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the  analysis  appears  accurate  within  the  bounds  of  existing  fracture 
mechanics  concepts. 

The  analysis  was  used  in  a  parametric  study  of  the  effects  of 
boron/epoxy  composite  reinforcement  on  crack  propagation  in  aluminum 
sheets.  The  study  showed  that  the  aspect  ratio  of  the  debond  area 
has  a  significant  effect  on  the  crack  propagation  in  the  aluminum 
sheet.  For  small  debonds  the  crack  propagation  rate  is  reduced 
significantly,  but  these  small  debonds  have  a  strong  tendency  to 
enlarge.  Debond  growth  is  most  likely  to  occur  in  reinforced  systems 
that  have  a  cracked  metal  sheet  that  is  reinforced  with  a  relatively 
thin  composite  sheet. 

The  analysis  can  be  used  to  predict  crack  groivth  in  reinforced 
systems.  Hence,  the  analysis  can  be  applied  in  developing  methods 
to  repair  damaged  metal  structure  and  to  increase  lives  and  payloads 
of  metal  structures  by  selective  reinforcement. 


APPENDIX  A 


DETERMINATION  OF  DEBOND  CONSTANTS 

As  discussed  in  Chapter  II,  debonding  can  be  predicted  in  a 
bonded  system  with  an  equation  of  the  type 

nj 

db/dN  =  C2  (G)  (5) 

where  G  is  the  strain  energy  release  rate  and  and  n2  are 

empirical  constants.  The  objective  of  this  appendix  is  to  determine 
the  empirical  constants  for  the  reinforced  system  used  in  the 
experimental  portion  of  this  dissertation. 

Specimen  Fabrication 

To  determine  Cj  and  several  test  specimens  with  the 

configuration  shown  in  figure  A.l  were  fabricated.  The  specimens 
consisted  of  1-inch  wide  strips  of  0.188  inch  thick  2024-T3  aluminum 
bonded  to  0.03  inch  thick  unidirectional  boron/epoxy.  The  strips  were 
bonded  with  Shell  EA-934  room  curing  adhesive.  To  maintain  a  constant 
adhesive  thickness  in  the  bond,  2  percent  by  volume  of  0.004  inch 
diameter  glass  beads  were  added  to  the  adhesive  prior  to  bonding. 

The  process  used  to  bond  the  aluminum  to  the  composite  was  as 
follows 
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Debond  specimen  configuration 
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SURFACE  PREPARATION 
Aluminum  2024-T3 

1.  Vapor  degrease  -  perch loroethylene  condensing  vapors  for 

5  to  10  minutes. 

2.  Grit  blast  with  220  grit  at  90  psig. 

3.  Alkaline  degrease  -  Oaklite  164  solution  (9  -  11  ounces/ 

gallon  of  water)  at  190  ±10°  F  for  15  minutes.  Rinse 
immediately  in  large  quantities  of  cold  running  water. 

4.  Acid  etch  -  place  panels  in  the  following  solution  for 

10  minutes  at  150°  ±5°  F. 

Distilled  water  30  parts 

Sulfuric  acid  (cone)  10  parts 
Sodium  Di chromate  1  part 

5.  Rinse  -  rinse  panels  in  clear,  deionized  running  water. 

6.  Dry  -  air  dry  15  minutes;  force  dry  10  minutes  at  150°  F 

±10°  F. 

Boron/epoxy 

1.  Vapor  degrease  as  above. 

2.  Grit  blast  with  220  grit  at  30  psig. 

BONDING 

1.  Bond  within  4  hours  of  surface  preparation. 

2.  Coat  surfaces  of  both  adherends  prior  to  bonding. 

3.  Cure  at  room  temperature  under  15  psig  ±2  psig  pressure. 

4.  Record  date  of  bonding. 
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The  method  of  surface  preparation  for  the  aluminum  was  taken  from 
Cagle  (1973)  while  the  surface  preparation  for  the  boron/ epoxy  and  the 
bonding  method  was  developed  by  the  author.  The  bonding  process  was 
verified  with  lap-shear  strength  tests  of  the  bonded  system. 

As  shown  in  figure  A.l  photoelastic  material  was  bonded  to  the 
surface  of  the  boron/epoxy.  The  photoelastic  material  enabled  tracking 
of  the  debond  front  in  the  fatigue  tests  of  the  specimens. 

Fatigue  Tests 

The  fabricated  specimens  were  tested  in  a  servo-hydraulic  fatigue 
machine  with  a  maximum  load  capacity  of  10,000  pounds.  All  of  the 
tests  were  conducted  at  a  loading  frequency  of  10  Hz  with  a  ratio  of 
minimum  to  maximum  load  in  the  load  cycle  of  R  =  0.05.  Duplicated 
tests  were  conducted  for  maximum  loads  of  5,000,  4,000,  and  3,000 
pounds . 

As  the  specimens  were  tested,  a  debond  developed  at  the  change 
of  cross  section  and  propagated  between  the  aluminum  and  boron/epoxy 
adherends.  Throughout  the  tests  the  location  of  the  debond  front 
was  indicated  by  an  isochromatic  that  was  observed  by  viewing  the 
photoelastic  material  through  a  polarizing  material.  The  location 
of  the  debond  front  is  plotted  against  the  number  of  applied  load 
cycles  on  figure  A. 2  for  all  of  the  tests.  The  results  of  these 
tests  will  be  used  with  a  stress  analysis  to  determine  the  empirical 

constants  c^  and  nj. 


1 OOK  200K 

applied  load  cycles,  N 
A. 2.  Debonding  behavior 
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Stress  Analysis 

Because  of  the  change  of  cross  section  of  the  test  specimen, 
as  an  axial  load,  P,  is  applied  to  it,  it  bends.  As  shown  by 
Timoshenko  (1961)  the  equilibrium  equation  for  a  beam  that  exhibits 
both  axial  and  bending  deformation  can  be  written  as 

dM  dw 

V  =  —  +  P —  (A.l ) 

dy  dy 

where  V  is  the  shear,  M  is  the  bending  moment,  and  P  is  the 
tensile  axial  load.  To  use  equation  (A.l)  both  the  moment  M  and 
the  axial  load  P  were  related  to  the  deflection  in  the  z-direction,  w, 
with  two  equations  given  by  Calcote  (1969)  as 


p  dv  d^w 

-  =  A  —  -  B,, - 

A  ”dy  '  dy^ 


dv  d^w 


M  =  B 


(A. 2) 


(A, 3) 


where  A  is  the  cross  sectional  area  of  the  beam  and  dv/dy  is  the 
axial  strain  of  the  beam  and 
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k  ^*^22 

An  =  ^  =-  E  ^  z'_i) 

i=1  ^  i=l 


1  k 


=:  Qll^i  -  2i-i>  '>k  '77 
•  ^  * 


(Ey)i 

•  (v  V  ) . 

^  xy  yx'^ 


where  i  indicates  the  layer  of  a  beam  with  k  layers  and  z  is 
measured  from  any  reference  surface. 

Eliminating  dv/dy  from  equations  (A. 2)  and  (A. 3)  yielded  an 
expression  for  the  moment  as 


d^w] 

d^w 

— 

P  +  B..  - > 

-  D  - 

A  ) 
11  ( 

11  dy2 

(A.4) 


Substituting  equation  (A.4)  into  equation  (A.l)  and  differentiating 
once  with  respect  to  y,  yielded  a  governing  equation  as 

(Bjj  1  d'*w  d^w 

Q  =i -  -D,,i -  +P  -  (A-5) 


-  Dji  > -  +  P  - 

( dy**  dy^ 


where  q^  is  a  transverse  distributed  load  acting  on  the  beam. 

The  boundary  conditions  for  equation  (A. 5)  were  determined 
from  the  end  conditions  of  the  test  specimen  installed  in  the  test 
machine  as  shown  schematically  on  figure  A. 3a.  With  the  assumption 


105 


that  the  specimen  is  effectively  fixed  at  both  ends  by  the  relatively 
massive  test  machine  grips  the  boundary  conditions  were  determined  as 


w(-d)  =  0  w(L)  =  0 

(A. 6) 

dw(-d)  dw(L) 

- =  0  -  =  0 

dy  dy 

At  the  change  in  cross  section  of  the  test  specimen  a  local 
moment  given  as 


M  =  Pe  {A-7) 

0 

where  e  is  the  distance  between  centroids  of  the  two  cross  sections, 
0 

was  produced.  With  this  moment  the  beam  can  be  modeled  as  a  symmetric 
beam  with  two  different  cross  sections  acted  upon  by  a  local  moment 
M  at  the  change  in  cross  section  and  an  axial  load  P  as  shown  in 
figure  A. 3b.  For  the  model  shown  in  figure  A. 3b,  equation  (A. 5) 
was  solved  with  conventional  finite  difference  techniques  (Ames  1971). 

To  verify  the  analysis,  the  strains  were  determined  by  both  the 
finite  difference  analysis  and  by  experiments.  The  test  specimen 
shown  in  figure  A. 3a  was  analyzed.  In  cross  section  2  the  specimen 
was  composed  of  four  layers;  the  metal  core,  the  adhesive  layer,  the 
composite  cover,  and  the  photoelastic  material.  The  thickness,  moduli, 
and  Poison's  ratio  for  each  of  these  layers  are  given  as 
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layer 

material 

thickness 

(in) 

elastic 

modulus 

psi 

Poisson's 

ratio 

1 

aluminum 

0.188 

10^ 

0.30 

2 

adhesive 

0.004 

600,000 

0.40 

3 

composite 

0.030 

3  X  10^ 

0.21 

4 

photoelastic 

0.083 

390,000 

0.36 

With  the  preceding  values  and  the  surface  of  the  metal  with  strain 
gages  (see  figure  A. 4)  as  a  reference  plane  B^,  and 

were  found  for  section  1  (see  figure  A. 3)  as 

A  =  2.1  X  10^  lbs/ in  B  =  200,000  lbs 
11 

D  =  2,500  lbs/ in 
1 1 

and  for  section  2  as 

Ajj  =  3  X  10®  Ibs/in  B^^  =  392,000  lbs 

Djj  =  65,000  Ibs/in 

With  the  use  of  the  previous  values  and  an  applied  load  P  of  5,000 
pounds  as  shown  on  figure  A. 4,  the  finite  difference  method  was  used 
to  calculate  deflections  and  curvatures  of  the  beam.  With  the  curva¬ 
tures  the  strains  in  the  beam  were  calculated  on  the  surface  of  the 
metal  with  (Calcote  1969) 

d^w) 


composite  ^ —  photoelastic 
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Comparison  of  finite  difference  solution  with  experimental  data 
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On  figure  A. 4  the  solid  line  indicates  calculated  strains  on  the 
surface  of  the  metal  adherend  while  the  circles  indicate  strains 
obtained  experimentally  with  strain  gages.  As  evident  from  the  figure, 
the  comparison  is  very  good.  Consequently,  the  developed  stress 
analysis  is  adequate  and  can  be  used  to  calculate  the  strain  energy 
release  rates  as  the  test  specimens  debond. 

Calculation  of  Strain  Energy  Release  Rates 
The  strain  energy  release  rate  can  be  calculated  as 


G 


lim 

\  AW 

AU 

Ab  ->•  0 

)  Ab 

Ab 

(A.  8) 


where  W  is  the  change  in  work  done  on  the  system,  U  is  the  change 
in  internal  strain  energy,  and  Ab  is  a  small  extension  of  the 
debond.  With  the  assumption  that  the  debond  extends  at  the  maximum 
applied  load,  the  work  done  as  the  debond  extends  can  be  calculated  as 


i’ 

AW  =  pjA6 


axial 


+  A6 


bending 


1 


+  MA00 


(A. 9) 


where  A00  is  the  rotation  of  the  beam  at  the  debond  front.  The 
change  in  internal  energy  can  be  calculated  as 


”  '^^axial  ^^bending 


(A. 10) 
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Substituting  equations (A. 9)  and  (A. 10)  into  equation  (A. 8)  and 
rearranging  slightly  yields  the  strain  energy  release  rate  as 


G  =  lim  <P 
Ab  ^  0 


^axial 


AU  .  , 
axial 


(A. 11) 


The  first  two  terms, which  neglect  bending,  of  equation  (A.ll)  can  be 
calculated  as  shown  by  Roderick  et  al  (1975)  as^  (for  a  unit  width) 


lim 

Ab  ->  0 


axial 


axial 


t  E 
c  c 


2E_r  (t,,,E  +  t-E  ) 
m  tn  111  111  c  c 


(A. 12) 


The  last  three  items,  which  are  due  to  bending,  can  be  calculated  with 
the  finite  difference  results  in  the  following  manner. 

First,  the  axial  deflection  due  to  bending  is  calculated  as 
shown  by  Den  Hartog  (1952)  as 


(A. 13) 


The  flexible  adhesive  and  photoelastic  material  are  neglected 
because  they  have  little  effect  on  the  strain  energy  due  to  the  axial 
load. 


A 


no 


wherG  thG  slope  dw/dy  was  calculated  at  each  nodal  point  in  the 
finite  difference  approximation.  Between  the  nodes  the  slope 
was  assumed  to  vary  linearly.  The  integration  of  equation  (A. 13) 
was  done  piecewise  over  the  length  of  the  beam.  The  work  done 
by  the  moment  was  expressed  as 

M  C  dw") 

H  .  -  (A-l'l) 

2  [dyj 

The  strain  energy  due  to  bending  was  calculated  as 

- dy  (A. 15) 

where  again  the  integration  was  performed  in  a  piecewise  fashion. 

Equations  (A. 13)  through  (A. 15)  were  evaluated  before  and  after 
a  debond  extension  Ab.  The  quantities  were  subtracted  to  calculate 
the  last  terms  of  equation  (A. 11).  The  result  was  added  to  equation 
(A. 12)  to  give  the  strain  energy  release  at  the  debond  front.  In 
figure  A. 5  the  solid  lines  show  the  calculated  strain  energy  release 
rate  plottedagainst the  debond  lengths  for  values  of  applied  loads  of 
5,000,  4,000,  and  3,000  pounds.  The  dashed  line  on  the  figure  shows  the 
strain  energy  release  rate--neglecting  bending--calculated  with 
equation  (A. 12)  As  evident  from  the  figure,  when  the  debond  length 
is  greater  than  0.5  inch  or  less  than  4.5  inches,  the  contribution 
of  bending  to  the  strain  energy  release  rate  is  small.  Particularly, 
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=  /' 


strain  energy  release  rate,  in-lb/in 


111 


4.0 


3.0 


total  strain  energy  release  rate 


- axial 


strain  energy  release  rate 


1 

1 


- 1 - 1 _ L 

2  3  4 

debond  length,  inches 


J. 

5 


Fig.  A. 5.  Strain  energy  release  rate  versus  debond  front 
location 


112 


when  the  debond  length  is  about  2,5  inches  long,  the  bending  contribution 
is  zero.^  Consequently,  for  a  debond  length  of  2. 5  inches, equation 
(A. 12)  can  be  used  to  calculate  the  strain  energy  release  rate. 

Curve  Fit  for  Empirical  Constants 

The  debond  rates  for  the  configuration  analyzed  were  determined 
by  taking  the  slope  of  the  curves  at  a  2.5-inch  debond  length  from 
figure  A. 2.  These  experimental  rates  and  the  corresponding  values 
of  the  calculated  strain  energy  release  rates  are  shown  in  Table  A.l. 

The  values  in  Table  A.l  were  used  to  determine  the  empirical  constants 
Cj  and  n^  in  equation  (5)  with  a  least  squares  fit  (Wylie  1966). 

To  perform  the  fit,  equation  (5)  was  written  as 


log(db/dN)  =  logCc^)  +  n2log(G) 


(A. 16) 


As  a  result  of  the  curve  fit,  Cz  and  nz  were  found  as 


c  =  3.158  X  10“^  n  =  3.616 

2  2 

Figure  A. 6  shows  the  data  as  points  and  the  fitted  equation  (5) 
as  a  solid  line.  As  evident  from  the  figure,  the  equation  fits  the 


^Actually,  there  are  also  two  other  debond  lengths  where  the 
strain  energy  release  rate  is  zero,  but  they  are  located  closer  to 
the  ends  of  the  specimen  where  the  analysis  may  be  more  inaccurate 
than  at  the  2.5-inch  debond  length. 
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TABLE  A.l 

STRAIN  ENERGY  RELEASE  RATES 
AND  DEBOND  PROPAGATION  RATES 


load 

lbs 

strain  energy 
release  rate,  G 
in-1 bs/in 

debond  rate 
db/dN 
in/cycle 

5,000 

2.15 

5.60  X  10'^ 

5,000 

2.15 

5.20  X  10"^ 

4,000 

1.38 

1.04  X  10"^ 

4,000 

1.38 

8.00  X  10"^ 

3,000 

0.77 

1.26  X  10"^ 

3,000 

0.77 

1.33  X  10"^ 

debonding  rate,  in/cycl 
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Fig.  A. 6.  Correlation  of  debonding  rates  with  strain  energy 
release  rates 
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data  well.  With  the  determined  empirical  constants,  equation  (5) 
can  be  used  to  predict  debond  propagation  rates  whenever  the  strain 
energy  release  rate  can  be  determined. 


APPENDIX  B 


ADHESIVE  SHEAR  DEFORMATION  ASSUMPTION 

As  mentioned  in  Chapter  IV,  the  complexity  of  the  analysis 
for  a  reinforced  system  is  significantly  reduced  by  assuming  that  the 
adhesive  only  undergoes  shear  deformation  while  the  adherends  only 
undergo  extensional  deformations  that  do  not  vary  through  the 
adherend  thickness.  Herein,  the  validity  of  these  assumptions  were 
investigated  by  comparing  a  one-dimensional  solution  that  was  based 
on  these  assumptions  with  a  more  rigorous  two-dimensional  finite 
element  solution.  Before  the  one-dimensional  and  two-dimensional 
solutions  were  compared  the  bulk  properties  of  the  adhesive  were 
determined. 


Adhesive  Bulk  Properties 

To  determine  the  bulk  properties  of  the  adhesive  an  appropriate 
test  specimen  was  designed  and  fabricated  in  several  steps.  First, 
a  female  plastic  mold  was  made  from  the  specimen  shown  in  figure  B.la. 
Then,  the  adhesive  liquid  base  and  hardener  were  combined  and  cast 
into  the  mold.  Next,  after  curing  24  hours  in  the  mold  the  adhesive 
specimen  was  removed  and  cured  an  additional  5  days  before  it  was 
handled.  Finally,  x-rays  were  taken  of  the  specimen  to  locate  air 
bubbles  developed  in  the  molding  process.  Specimens  that  contained 
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large  air  bubbles  in  critical  areas  were  scrapped.  Six  test 
specimens  were  fabricated  in  this  manner  but  only  three  were  acceptable 
for  testing. 

Prior  to  the  testing  the  specimens  were  instrumented  with  strain 
gages  and  linear  variable  differential  transformers  (LVDT)  as  shown 
in  figure  B.lb.  The  strain  gages  were  used  to  obtain  the  longitudinal 
and  transverse  strain  while  the  LVDT's  were  used  to  check  the  strain 
gages  (the  LVDT's  were  wired  to  eliminate  bending  effects  in  their 
readings  while  the  strain  gages  would  include  bending  in  their 
readings)-  Discrepancies  between  the  readings  would  indicate  bending 
due  to  poor  specimen  alignment  in  the  test  machine.  Each  instrumented 
specimen  was  placed  in  a  servo- hydraul ic  test  machine  with  a  loading 
range  of  2,000  pounds  and  a  sensitivity  of  ±10  pounds.  Then,  each 
specimen  was  loaded  to  failure  at  a  rate  of  80  Ibs/sec. 

In  this  manner,  the  three  specimens  were  tested  and  the  results 
were  nearly  identical  for  all  three  of  the  tests.  For  all  specimens, 
the  longitudinal  strain  calculated  from  LVDT's  agreed  within  1.5  percent 
of  the  longitudinal  strain  gage  reading  and  indicated  the  specimens 
were  aligned  properly  in  the  test  machine.  The  data  obtained  from  the 
strain  gages  are  shown  on  figure  B.2  in  the  form  of  a  stress-strain 
plot.  On  the  figure  the  heavy  solid  line  indicates  the  strain  in  the 
loading  direction  Cy  flnd  the  dashed  line  indicates  the  strain 
transverse  to  the  loading  direction  e  .  As  indicated  by  the  light 
solid  lines  on  the  figure,  the  stress-strain  curve  can  be  approximated 
by  a  bilinear  stress-strain  curve  with  a  change  of  slope  occurring 


at  4,200  psi. 


stress,  ksi 


Fig.  B.2.  Adhesive  stress-strain  curve 
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In  the  initial  linear  region  the  following  values  were  obtained 
from  the  curve: 

Aa=  3,000  psi,  Ae=  0.005,  e  =  0.002 

y  y 

Using  these  parameters,  the  Poisson's  ratio  and  the  elastic  modulus 
were  calculated  as 


V  =  —  =  0.40  E  =  - —  =  600,000  psi 

e  y  Ac, 

y  y 

With  the  assumption  that  the  adhesive  is  isotropic  the  shear  modulus 
was  calculated  as 


e  =  -  =  215,000  psi 

2  (1  +  V  ) 

yx' 

Along  similar  lines,  the  material  parameters  in  the  second  linear 
region  were  determined  as 

V  =  0.28  E  =  190,000  psi  G  =  74,000  psi 

xy  y 

The  fracture  of  the  specimen  occurred  at  a  stress  of  6,600 
psi  at  an  axial  strain,  e^,  of  0.0215.  With  the  preceding  material 
property  values  for  the  adhesive,  the  one-dimensional  and  two- 
dimensional  solutions  were  compared. 
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As  a  test  case,  the  shearing  stresses  in  the  adhesive  layer 
of  the  specimen  shown  on  figure  B.3a  (the  specimen  is  symmetric  about 
the  x-y  plane)  was  calculated  with  both  types  of  solutions  for  a  plane 
stress  state.  Although  the  solution  was  for  plane  stress,  for  the 
self-equilibrating  load  system  shown  on  figure  A. la  the  stress 
distributions  are  identical  for  both  plane  stress  and  strain  states 
(Timoshenko  (1951)).  Consequently,  although  the  tests  case  considered 
a  state  of  plane  stress,  the  results  are  also  applicable  to  a  state  of 
plane  strain  which  may  be  more  appropriate  for  a  section  taken  through 
the  thickness  of  the  reinforced  system  shown  on  figure  9. 

One-Dimensional  Solution 

Figure  B.3b  shows  a  freebody  of  the  specimen  shown  in  figure  B.3a. 
On  the  figure  P  is  half  the  load  in  a  composite  adherend,  F(y) 
and  q(y)  are  the  load  in  the  metal  adherend  and  the  shear  flow  in 
the  adhesive  at  any  point  y,  and  t^^,  2t^,  and  t^^  are  the 

thicknesses  of  the  metal,  composite,  and  adhesive  respectively. 

The  change  in  the  load  F(y)  with  respect  to  y  is  the  shear 
flow  in  the  adhesive  layer  given  as 

dF(y) 

q(y)  =  - 

dy 

The  shear  stress  in  the  adhesive  is  this  shear  flow  divided  by  the 
width,  w,  of  the  specimen,  given  as 

q  1  dF(y) 

T(y)  =  —  = -  (B.2) 

w  w  dy 


adhesive 


L 


composite 


(b)  Freebody  for  one-dimensional  solution 


Fig.  B.3.  Specimen  configuration  and  freebody  for 
one-dimensional  solution 


With  the  assumption  that  the  adhesive  only  undergoes  shear  deformation 
and  that  this  deformation  does  not  change  through  the  adhesive  thickness, 
the  shearing  stress  in  the  adhesive  can  be  related  to  the  shearing 
strain  in  the  adhesive  by  the  constitutive  equation  as 

T(y)  =  Y(y)Gg,j  (B.3) 

The  shearing  strain  in  the  adhesive  can  be  related  to  the  extensional 
displacements  in  the  metal,  composite,  >^Q(y)» 

as 

V  (y)  -  V  (y) 

Y{y)  = -  (B-4) 

^ad 

Substituting  equation  (B.4)  into  equation  (B.3)  and  that  result  into 
equation  (B.2)  yields 

dF(y)  GgjW 

-  =  -  fv,(y)  -  v^(y)]  (B.5) 

‘‘J'  ‘ad 


Equation  (B.5)  can  be  differentiated  with  respect  to  y  to  yield 
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when  the  derivatives  of  the  extensional  displacements  and 
with  respect  to  y  are  the  extensional  strains  in  the  metal  and 
composite  adherends  respectively.  These  strains  were  related  to  the 
extensional  loads  in  the  adherends  by 


dv„(y) 

F(y) 

dv^(y)  P  -  F(y) 

dy 

"'m'm 

(B.7) 


where  E|^  and  are  the  moduli  of  the  metal  and  composite  respec¬ 
tively.  Substituting  equations  (B.7)  into  equation  (B.6)  yielded  a 
second  order,  nonhomogeneous  differential  equation  in  F(y)  as 

d^F(y) 

-  -  aF(y)  =  3  (B.8) 

dy 

where 


where  s  is  the  stress  in  the  composite  and  is  the  shear  modulus 

of  the  composite. 

Solving  equation  (B.8)  yielded  a  complete  solution  as 
Ay  -Ay  3 

F(y)  =  Cje  +026  +  -  (B.9) 

a 

With  referral  to  figure  B.3a,  the  boundary  conditions  for  equation  (B.9) 
are 


y  =  0  F(y)  =  F(0)  =  0 

y  =  00  F(y)  is  bounded 


(B.IO) 


With  the  use  of  these  boundary  conditions,  the  constants  in  equation 
(B.9)  were  found  to  be 

3 


By  the  use  of  equation  (B.2),  the  shear  stress  in  the  adhesive  was 
calculated  from  equation  (B.ll)  as 
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T(y)  = 


t  F 


e 


-/ay 


(B.12) 


Finite  Element  Solution 

A  finite  element  computer  program,  PLANE,  was  used  to  calculate 
the  shearing  stresses  in  the  adhesive  layer  in  the  specimen  shown  in 
figure  B3.a.  PLANE  is  an  elastic-plastic,  two-dimensional  program 
which  uses  constant  strain  and  linear  strain  triangular  elements. 

PLANE  was  developed  by  the  Grumman  Aircraft  Corporation  for  the 
National  Aeronautics  and  Space  Administration  and  was  documented  by 
Armen  and  Levy  (1962).  The  mesh  used  for  the  analysis  is  shown  on 
figure  B.4  and  contains  1,522  degrees  of  freedom.  The  triangles 
are  predominately  linear  strain  triangles  which  allow  linear  variations 
in  the  stresses  and  strains  through  the  elements.  Each  of  the  adherends 
is  modeled  with  several  elements  through  the  thickness  thus  allowing 
variations  of  extensional  and  shearing  stresses  through  the  thickness. 

In  contrast,  the  one-dimensional  analysis  assumed  uniform  extensional 
stresses  and  no  shearing  stresses  through  the  thickness  of  each 
adherend.  The  adhesive  layer  was  modeled  by  one  layer  of  elements 
which  allowed  linear  variation  in  stresses  through  the  adhesive 


thickness. 
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element  mes 
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One-Dimensional  Versus  Finite  Element  Solution 

To  compare  the  results  of  the  two  solutions,  the  specimen 
configuration  shown  in  figure  B.3a  with  a  width  of  w  =  1.0  inch 
and  the  following  parameters 

metal  adhesive  composite 

thickness  t^  =  0.1  in  t^^  =  0.004  in  t^  =  0.1  in 

moduli  =  10.3  X 10®  psi  Gjj  =  215,000  psi  E^  =  30  x  10®  psi 

and  an  applied  stress  s  of  1,000  psi  was  analyzed  with  both  solutions. 
The  parameters  used  in  these  solution  were  typical  for  constituents 
used  in  the  reinforced  system  discussed  in  Chapter  II.  The  adhesive 
shear  stress  calculated  from  the  one-dimensional  solution  equation 
(B.12)  was  plotted  as  a  solid  line  against  y  (distance  from  the 
edge  of  the  metal  adherend)  on  figure  B.5.  On  the  same  figure  the 
circles  indicate  values  of  the  adhesive  shear  stress  calculated 
from  the  finite  element  solution.  As  evident  from  the  figure  the 
one-dimensional  solution  gives  shear  stress  values  twice  as  high  as 
those  obtained  from  the  finite  element  solution  near  the  edge  of  the 
metal  adherend  (y  =  0).  Evidently,  the  shearing  deformation  of  the 
adherends  which  was  not  accounted  for  in  the  one-dimensional  solution 
has  a  significant  effect  on  the  values  of  the  shear  stresses. 

To  account  for  the  adherend  shear  deformation  and  still  use  the 
simplified  one-dimensional  analysis,  an  effective  shear  modulus 


shear-stress  in  adhesive 
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Fig.  B.5.  Finite  element  versus  one-dimensional  solution 
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was  introduced.  The  magnitude  of  was  determined  by  equating 

the  maximum  shear  stresses  in  the  adhesive  calculated  by  the  finite 
element  solution  to  the  expression  for  the  maximum  shear  stress  from 
the  one-dimensional  solution  from  equation  (B.12)  as 


T 

max 

Finite  Element 


(B.13) 


Solving  equation  (B.13)  for  6  yielded  an  expression  for  an 
effective  shear  modulus  as 


With  this  value  for  G  in  equation  (B.12)  yields  a  corrected  one¬ 
dimensional  solution  for  the  adhesive  shear  stress  as 

T(y)  =  - =— 

‘adEc 

For  the  sample  analysis,  equation  (B.14)  yielded  Gg^^  as  64,000 
psi.  Using  this  value  for  Ggff>  equation  (B.15)  was  plotted  against 
y  as  a  dashed  line  on  figure  B.5.  The  agreement  between  equation  (B.15) 
and  finite  element  results  indicated  by  the  circles  on  the  figure  is 
excellent.  Consequently,  the  assumptions  made  in  Chapter  IV,  that  the 
adherends  only  undergo  extensional  deformation  while  the  adhesive 
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only  undergoes  shear  deformation,  can  be  used  to  accurately  predict 
shear  stresses  in  the  adhesive  with  a  simplified  analysis  if  an 
effective  shear  modulus  for  the  adhesive  is  used  in  the  calculations. 

For  analysis  of  the  reinforced  system  shown  in  figure  2  the 
effective  modulus  was  determined  for  a  range  of  adherend  thicknesses 
for  the  adhesive  used  in  the  reinforced  system  both  for  before 
and  after  the  adhesive  yields.  To  determine  the  values  of  effective 
shear  moduli  numerous  finite  element  solutions  were  run  with  different 
adherend  thicknesses  for  both  the  initial  bulk  adhesive  shear  modulus 
of  215,000  psi  and  the  bulk  shear  modulus  after  yielding  of  74,000 
psi,  The  results  of  these  calculations  are  given  in  table  B.l.  In 
the  table  the  maximum  shear  stress  calculated  from  the  finite  element 
solution  and  the  effective  shear  modulus  are  tabulated  for  the 
different  adhesive  thicknesses. 

As  shown  in  table  B.l  the  value  of  the  effective  modulus  for  the 
initial  shear  modulus  of  the  adhesive  does  not  vary  much  with  adherend 
thicknesses.  The  average  value  for  6^^^  is  about  65,000  psi  and 
is  within  3  percent  of  any  of  the  calculated  values. 

Also,  as  shown  in  table  B.l,  the  value  of  the  effective  modulus 
for  shear  modulus  of  the  adhesive  after  yielding  also  has  little 
variation  with  adherend  thicknesses.  The  average  value  in  this  case 
is  about  36,000  psi  and  is  within  3.3  percent  of  any  of  the  calculated 
values. 

As  evident  from  the  previous  discussion  the  reinforced  system 
can  be  analyzed  by  assuming  that  the  adherends  only  undergo  extensional 
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deformation  while  the  adhesive  only  undergoes  shearing  deformation 
if  an  effective  shear  modulus  of  65,000  psi  is  used  for  the  adhesive 
before  the  adhesive  yields  and  an  effective  shear  modulus  of  36,000  psi 
is  used  for  the  adhesive  after  yielding. 


EFFECTIVE  SHEAR  MODULI  FOR  REINFORCED  SYSTEM 
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'  ,  psi  671  876  1,002  1,090 

fTla  X 

psi  36,840  36,000  35,400  35,000 


APPENDIX  C 


REMOTE  STRESSES  IN  THE  ADHERENDS 

For  equilibrium  to  exist,  the  macroscopic  stresses  applied  to  the 
reinforced  system  must  be  balanced  by  the  stresses  in  the  adherends 
of  the  system.  This  relationship  can  be  written  in  vector  form  as 


(C.l) 


where  A  represents  the  area  of  the  reinforced  system  and  A^  and 
A^  represent  the  area  of  each  element.  For  the  case  where  stress 
was  applied  in  only  the  y-direction, 


With  strains  uniform  through  the  thickness,  (C.l)  can  be  rewritten  for 
a  unit  width  as 
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0 

V 

(tm  t_)  = 

'  m  c 

[C]t„  +  [D]t, 

"y 

0 

S. , 

xy 

{C.3) 


where  the  C  and  D  are  the  stiffness  matrices  for  an  isotropic 
and  orthotropic  material  in  plane  stress,  respectively,  as  given  by 
Zienkiewicz  (1971)  as 

E  vE 

-  -  0 

1  -  1  - 


[C] 


vE  E 

-  -  0 

1  -  1  - 


(C.4) 


[D]  = 


V  E 
yx  X 


^yx^x 


^yx 


(1 


2 

y 


1  - 


yx 


(C.5) 
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with  E,  G  ,  and  V  „  denoting  the  extensional  modulus,  the 
yx  yx 

shear  modulus,  and  Poisson's  coefficient  respectively.  With  simple 
matrix  algebra,  equation  (C.3)  can  be  solved  and  the  strains  determined. 
These  strains  in  turn  can  be  used  with  the  appropriated  stiffness 
matrices  to  calculate  the  stresses  in  the  metal  and  composite 
adherends  as 


■^^x 

=  [C] 

^"’yy 

{C.6) 

;^y 

,"^y 

’'^Sx 

’"Sx' 

ac 

yy 

=  [D] 

(C.7) 

0C 

xy 

."V 

APPENDIX  D 


GREEN'S  FUNCTION  FOR  THE  CRACKED  SHEET 

The  solutions  for  the  interlaminar  stresses  and  the  strain  energy 
release  at  the  debond  front  developed  in  Chapters  IV  and  V  respectively 
require  Green's  functions  for  both  displacements  and  stresses.  As 
pointed  out  by  Dennemeyer  (1968),  the  solution  for  a  concentrated 
load  acting  on  a  body  can  be  used  as  a  Green's  function. 

Herein,  the  solution  for  four  concentrated  loads  that  are 
symmetric  with  respect  to  a  crack  in  an  isotropic  sheet  was  developed. 

The  solution  was  based  on  elasticity  theory  using  complex  variable 

theory  as  developed  by  Muskhelishvili  (1975).  The  solution  was 
predicted  on  the  assumptions  that  the  cracked  sheet  is  infinite 
isotropic,  and  can  be  described  by  a  plane  stress  or  strain  analysis. 

As  shown  by  Muskhelishvili  (1975),  both  the  stresses  and  displacements 
in  a  cracked  sheet  can  be  expressed  in  terms  of  two  stress  functions, 
$(z)  and  'i';z),  as 

=  Real{$(z)  +  «I>(z)  -  Cz$'(z)  +  'i'(z)]}  (D.l) 

Oy  =  Real{$(z)  +  1>(z)  +  z$' (z)  +  'f(z)}  (D.2) 
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a  =  Imag{z'l>' (z)  +  ’l'(z)}  (D-3) 

xy 

2G(u  +  iv)  =  n(!)(z)  -  (z)  -  »i^(z)  (D-^) 

where  a  ,  a  and  a  are  stresses,  u  and  v  are  displacements, 

X  y  xy 

the  bar  over  the  functions  denotesthe  complex  conjugate,  i  is  the 
square  root  of  -1  and 

d(j)(z)  #(z) 

$(z)  =  (})'(z)  =  -  4'iz)  =  il;'(z)  =  - 

dz  dz 

=  3  -  4v  plane  strain 
3  -  V 

n  =  -  plane  stress 

1  +  V 

is  the  Poisson's  ratio. 

The  stress  functions,  ^(z)  and  '{'(z),  for  the  cracked  sheet 
under  four  concentrated  loads  as  shown  on  figure  D.lc  were  constructed 
by  superimposing  the  stress  functions  for  a  cracked  sheet  under  two 
different  loading  conditions.  The  first  condition  which  is  shown 
on  figure  D.la  has  concentrated  loads  acting  on  the  cracked  sheet  plus 
a  stress  distribution  applied  to  the  crack  surface.  This  stress 
distribution  was  equal  to  the  stress  distribution  which  exists  along 
the  x-axis  for  an  uncracked  sheet  with  concentrated  loads.  This  stress 
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ree 
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distribution  closes  the  crack  and,  therefore,  effectively  eliminates 
it.  As  a  result,  figure  D.la  represents  an  uncracked  sheet  with 
concentrated  loads  acting  on  it.  The  second  condition  which  is  shown 
on  figure  D.lb  is  for  a  stress  distribution  acting  on  the  surface  of  a 
crack.  This  stress  distribution  has  the  same  magnitude,  but  is 
of  opposite  sign  to  the  stress  distribtuion  applied  to  the  crack  in 
the  first  condition.  The  development  of  the  stress  functions  for 
these  two  conditions  follows. 

The  stress  functions  for  figure  D.la  were  developed  by  super¬ 
imposing  the  solution  for  single  concentrated  forces  that  act  in 
different  quadrants.  Fora  point,  z^,  in  the  first  quadrant  the 
solution  was  given  by  Muskhel ishvi  1 i  (1975)  as 


1 


0(z)  =  -S 


'i'(z)  =  Sn  - 
z 


(D.5) 


where 


X  +  iY 

S  =  - - — 

2tt(1  +  n)t|^ 

and  X  and  Y  are  the  load  components  in  the  x  and  y  directions 
respectively,  tt  =  3.1459,  and  t^^  is  the  thickness  of  the  sheet. 
Equations  (D.5)  were  generated  for  the  second,  third  and  fourth 
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quadrants  by  replacing  S  and  with  -S  and  -z^,  -S  and  -z^, 
and  ^  and  respectively.  These  stress  functions  were  then 
superimposed  to  form  the  functions 


(D.6) 


(0.7) 


where  upon  differentiation  equation  (D.6)  becomes 
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4»(z)  =  S  [Log(z  +  Zq)  -  Log(z  -  z^)  ] 
+  S  [Log(z  +  Zjj)  -  Log(z  -  z^)  ] 


(D.IO) 

Equations  (D.6)  through  (D.IO)  are  the  required  stress  functions, 
derivative,  and  integrals  to  compute  stresses  and  displacements  from 
equations  (D.l)  through  (D.14)  for  figure  D.la. 

The  stress  functions  for  load  condition  2  shown  on  figure  D.lb 
were  obtained  in  the  following  manner.  Following  Muskhelishvili  (1975), 
a  new  stress  function  fi(z)  was  introduced  which  is  related  to  the 
previously  discussed  stress  function  by  the  equation 

'l'(z)  =  n(z)  -  $(z)  -  z$'(z)  (D.ll) 

where 


q{z)  =  q{z) 


The  required  stress  functions,  $(z)  and  ^2(z),  can  be 
determined  for  pressure  acting  on  a  single  crack  as  (Muskhelishvili 


1975) 

1  I(t)p(t)dt  1  q(t)dt 

♦„(2)  - -  -  +—  -  (D.ia) 

P  2TriI(z)  -a  t  -  z  2TTi  -a  t  -  z 


1  I(t)P(t)dt  1  q(t)dt 

n(z)=— —  r - -  (0-13) 

P  2TTiI(z)  -a  t  -  z  2TTi  -a  t  -  z 


with 

crack  length 

I(z)  =  /z^  -  a^  a= - 

2 


a  and  T  are  normal  and  shearing  stresses  acting  on  the  crack 

y  xy 

surfaces  respectively  (the  plus  sign  indicates  the  upper  surface  of 
the  crack  while  the  minus  sign  indicates  the  lower  surface.)  P(t) 
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and  q{t)  are  total  normal  and  shearing  pressures  which  act  on  the 
crack  surfaces.  As  mentioned  previously,  the  stresses  which  acts 
on  the  crack  surface  in  figure  D.lb  is  equal  to  the  magnitude  but  of 
opposite  sign  to  the  stress  caused  by  four  point  loads  acting  on  an 
uncracked  sheet.  Therefore,  the  functions,  <i>(z)  and  v(z),  shown 
in  equations  (D.6)  and  (0.7),  which  are  solutions  for  point  loads  on  a 
continuous  sheet,  can  be  used  to  find  P{t)  and  q{t). 

The  normal  and  shearing  stresses  acting  on  the  crack  were 
calculated  (Muskhelishvili  1975)  in  terms  of  equations  (0.6),  (D.7),  and 
(D,8),  as 


a  -  iT  =  ^(z)  +  $(z)  +  z$'  (z)  +  T’(z)  (D.15) 

y  xy 

Substituting  equations  (0.6),  (0.7),  and  (0.8)  into  equation  (0.15) 
yields 


iT 


xy 


a{z,2^)  +  a{z,z^) 


(0.16) 


or 


Oy  -  iT^y  =  2Real  {a{z,z^)) 


(0.17) 


where 


z  +  z. 
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Equation  (D.16)  shows  that  no  imaginary  term  exists.  Consequently, 

T  is  zero.  Also,  a  has  the  same  value  along  the  x-axis  independent 

xy  y 

of  the  direction  in  which  the  x-axis  is  approached.  Consequently, 

0  '*’  =  a  Therefore,  with  consideration  given  to  the  above  statements 

y  y 

and  equations  (D.14),  the  total  pressures  on  the  crack  are 
P(t)  =  aita^)  +  a(t,ZQ) 

(D.19) 

q(t)  =  0 

As  a  result  of  the  above  simplifications,  equations  (0.12)  and  (D.13) 
become  equal  to  each  other  and  are  expressed  as 

1  I(t)P(t) 

${z)  =  f2(z)  =  -  /  ®  -  dt  (D.20) 

2iTiI(z)  -a  t  -  z 


The  integral  in  equation  (D.20)  was  evaluated  by  contour 
integration  along  the  contour  shown  in  figure  D.2  by  using  the 
Residue  theorem  given  as  (Wylie  1960) 
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Fig.  D.2.  Path  for  contour  integration 
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The  contour  shown  in  figure  D.2  can  be  broken  into  several  sections. 
Thus,  equation  (D.19)  becomes 


/  A(^)  +  /  A(c)  +  /  A(c)  -  /  A{?)  +  /  A(?) 

r,  r,  r, 


-  f  A(5)  +  f  A(c)  =  2T7i  E  Res  (D.23) 


where 

p(c) 

A(c)  =  -  dc,  C  =  t  +  iS 

C  -  z 


On  figure  D.2  the  integral  is  evaluated  along  the  contour  as  R 
approaches  infinity  and  e  approaches  zero.  Each  of  the  contours 
can  be  expressed  as  follows. 


Kc)  P(?) 

f  A(c)  =  lim  / -  dc 


e  ->  0 


Let  c  +  a  =  ce 


C  -  z 


d^  =  iee^®d0 


ie  Ve^e^®  -  2aee^®  P(ee^® 


a)e^®d0 


'•  /  A(c)  =  lim 

^1  e  -»•  0 


/  ‘ 
217 


-  a  -  z 


(D.24) 
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Noting  that 

+  - 

I(t)  =  I(t)  for  c  >  a 

the  sum  of  the  next  two  segments  of  the  contour  integral  equals 
zero,  because 

+  + 

I(t)  P(t)dt  I(t)  P(t)dt 

/  A(0  =  -  f  A(?)  =  - 

r  a  t-z  r  R  t-z 

5  6 

SO  that 

/  A(c)  +  /  A(0  =  0  (D.29) 

The  last  contour  segment  can  be  evaluated  as 

Kc)  P(c)d? 

/  A(c)  =  lim  /  - 

^  R  -  «  ^  " 

Let  c  =  Re^®,dc  =  Rie^®d0 

Ri  VR'e'''®  -  P(Re^®)d0 

•  2Tr 

••  /  A(c)  =  lim  /  - in - 

r,  „  .  _  •  -  z 
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Substituting  equations  (D.24),  (D.25),  (0.28),  (D.29),  and  (D.30) 
into  equation  (D.23)  and  rearranging  the  terms  yield  the  desired 


integral  as 


I(t)  P(t)dt 


2aee^®P(ee^® 


t  -  z 


=  Tri  E  Residues  -  lim  -  / 

r)  2Tr 

e  ->•  0 


ee^®  -  a  -  z 


1  ieVe^e^®  +  2aee^®  P(ee^®  -  a)e^®d0 

-TT 

-lim  -  /  - —Jq  “  ■ 

2  TT  ee'''  +  a  -  z 

£  ->  0 


1  RiUR^e^''®  -  a2  p(Re''=')d0 

21T  ' 

-  lim  -  /  - fs 

„  2  ^  Re  -  z 

R  ^  00 


(D.31) 


where  P(t)  is  given  by  equation  (D.19).  As  an  example,  equation 
(D.31)  was  calculated  for  the  first  term  of  P(t).  For  this  case 
the  expression  to  consider  was  chosen  as 


Kc)  P(C) 


K  -  z 


-4I(?)  Sz, 


-  Zq2)(^  -  z)  1 


,  P(t)  = 


(t^  -  ZQ2)(t  -  z) 


(D.32) 
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The  residuals  for  expression  (D.32)  were  calculated  with  equation  (0.22). 
They  become 


To  insure  single  valuedness  of  the  stress  functions,  values  of  1(c) 
must  be  chosen  so  that  they  lie  on  the  same  branch.  Values  of  1(c) 
will  lie  on  the  same  branch  for  all  values  for  c  if 


1(0 


=  1 


(D.34) 


The  two  possible  values  of  1(c)  are  the  complex  numbers  w^  and  -w^. 
Hence,  equation  (D.34)  requires  that  I(Zq)  equals  w^  and  that  I(-Zq) 
equals  -w^,  or  simply  that  Kz^)  =  -I(-Zq).  With  the  proper  values 
of  1(c)  defined  according  to  equation  (D.34),  equation  (D.33) 
reduces  to 
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Residual  =- 


"o  -  ^ 


-  (zKz^)  -  z^I(z)) 


(D.35) 


In  order  to  evaluate  equation  (D.31)  the  limits  of  three 
integrals  must  be  taken.  The  first  two  integrals  involve  limits  as 
e  0  (equations  (D.24)  and  (D.25)).  The  evaluation  of  these  limits 
is  similar  for  both  equations.  The  limits  can  be  obtained  by  rearranging 
the  terms  in  the  integrand  so  that  they  can  be  expressed  by  a  simple 
binomial  series  as 


(x  +  y)"^  =  x"^  +  nx'^”V  + 


n(n  -  1) 


<  x") 


(D.36) 


(1  ±  x)“^  =  1  +  X  +  x=^+  x^  +...  (x^  <  1) 


(D.37) 


With  the  use  of  the  first  term  of  P(t)  as  shown  in  (D.32),  as  an 
example  (D.24)  can  be  expressed 


-4zQie  -  2aee^®e^®d0 


/  A(?)  =  lim  S  /  — -g - -  ie” 

r  27T  (ee  ^  -  a  +  z  )(ee 

e  ->*  0 


a  -  ZQ)(ee  -  a  -  z) 


(D.38) 
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The  term  Vee^  -  2aee^®  in  equation  (D.39)  can  be  expressed  in  terms 
of  equation  (D.36)  as 


=  (1 

=  1  [, 


2a 


ee 


.i0 


2a 


r- 


ee 


ie 


4a 


e  e 


32a‘ 


(0.39) 


while  the  term  (ee^®  -  a  +  z)'^  can  be  expressed  in  terms  of  equation 
(0.37)  as 


(D.40) 
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Substituting  equations  (D.39),  (0.40)  and  similar  expressions  for  the 
remaining  terms  in  the  denominator  of  equation  (D.38)  into  equation 
(D,38)  yields  an  expression  as 


(D.41) 


whose  limit  is  zero.  Equation  (D.25)  can  be  evaluated  by  similar  means. 
The  result  for  P(t)  as  shown  by  (D.32)  is  also  zero. 

Equation  (D.30)  can  also  be  evaluated  by  a  limiting  process  and 
by  the  use  of  equations  (D.36)  and  (D.37)  as  R  For  this  case, 

the  integrand  has  to  be  arranged  in  a  manner  that  makes  the  expansions 
for  equation  (D.36)  and  {D.37)  valid  for  large  values  of  R.  For  exam¬ 
ple,  for  the  first  term  of  P(t)  as  given  by  (D.37),  equation  (D.30) 
can  be  expressed  as 
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“  Re^®l 


/  A(c) 

r 


(D.42) 


By  substituting  the  appropriate  expansions  for  the  terms  in  the  integrand 
of  equation  (D.42),  the  limit  as  R  increases  without  bound  gives 
the  value  of  equation  (D.42)  as  zero. 

Similarly,  equation  (D.31)  can  be  evaluated  for  all  tenns  of 
P(t)  contained  in  equation  (D.19).  All  of  these  terms  were  combined, 
simplified,  and  substituted  into  equation  (D.20)  to  yield 


$(z,Zq)  =  f2(z,Zj^)  =  S  B(z,Zq)  +  S  B(z,Zq) 


(D.43) 


where 


z  -  z. 


Z  +  Z, 


2  z^-z^^\^-v'(z-z^)^'(z.z^)^ 


^0  -  "o 


2I(zJI(z)  )  (z^  -  z)2  (z  +  z)2  ( 


KZq) 


Z  2  _  z2 


nl(Zo) 


(D.44) 


A 
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With  the  use  of  equation  (D. 43),  which  shows  4'(z)  and  n(z) 
to  be  equal,  and  equation  (D.ll),  H'^z)  was  expressed  as 

M^iz)  =  ${z)  -  $(z)  -  z4>'(z)  (D.45) 

where  <j)(z)  is  given  by  right  hand  side  of  equation  (0.43).  To 
evaluate  equation  (D.45),  the  derivative  of  $(z)  was  required.  Hence, 
the  derivative  of  B(z,Zq)  was  required  because 

d<I)(z) 

$'{z)  =  - =  SB'(z,Zq)  +  SB'(z,Zq)  (0.46) 

dz 


Using  a  symbolic  manipulation  system,  MACSYMA  (1975),  written 
in  LISP  programming  language,  the  derivative  of  B(z,Zq)  was  found 
to  be 


B'(z,Zq)  =  -2z 


(z^  - 


3^0%  * 


(z„  -  Z)>  (z„  +  z) 


1  f  1  rnl(?„)(a%“  -  2z'  +  a“z“) 

I(z)  jz"  -  a"  I  dj,"  - 

I(z„)(ad^“  -  2z'  a-  ad^) 

(z„"  - 

0  j 
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(D.47) 


To  complete  the  evaluation  of  the  functions  used  in  equations 
(D.l)  through  (0.4)  for  pressure  acting  on  the  crack  surface,  the 
integrals  of  $(z)  and  H'iz)  given  by  (l)(z)  and  ij;(z)  respectively 
were  determined.  As  shown  by  equation  (D.45),  once  (t)(z)  is  found 
ij;(z)  is  also  determined.  To  evaluate  (})(z)  by  integrating  (D.43) 
the  integral  of  B(z,Zq)  was  determined  (evaluation  of  the  integral 
of  B(z,z^)  is  the  same  as  for  the  integral  of  B(z,Zq)  except  that 
Zq  is  replaced  by 

Many  of  the  terms  in  B(z,Zq)  are  easy  to  integrate  by  using 
standard  rules  of  calculus.  However,  those  terms  which  contain  I(z) 
in  the  denominator  require  some  rearrangement  before  the  integration 
is  attempted.  For  example,  integrals  of  the  type 

dz 

l(z)(Zo  " 


Il(z)  =  / 


(D.48) 
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appear  frequently  and  were  integrated  as  follows.  First,  a  change 
of  variables  was  made  by  letting 

z  =  asinG,  dz  =  acosGdG 

With  the  restriction  that  |zj  <  a,  equation  (D.48)  became 

dz  acosGdG 

Il(z)  =  /  -  =  / - (D.49) 

I(z)(Zp  ~  z)  I(z)(Zq  -  asinG) 

Next,  the  trigonometric  relation  sin^G  +  cos^G  =  1  was  used 
to  express  cosg  as 

)I7~7 

cosG  =  -  (D.50) 

a 

At  this  point  care  must  be  taken  to  choose  the  correct  value  of  the 
multi-valued  function.  The  correct  value  was  assured  by  requiring 
equation  (D.34)  to  hold.  For  example,  I(z)  can  be  written  as 


I(z)  =  \/z^  -  a^  =  i  V  a^  -  z^  or  -  i\/a^  -  z^ 
but  according  to  equation  (D.34) 


1  im 


z2  -  a 


1  im 


-tiya 


=  1 


=  1 


2  oo 


Z 


Z  00 


Z 
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so  then 


I(z)  =  -  i  \/a^  -  or  \/a^  -  z^  =  il(z)  (D.51) 

Substituting  equation  (D.51)  into  equation  (D.40)  and  that  result 
into  equation  (D.49)  yields  the  result 

dz  d0 

f - =  i  / -  (D.52) 

I{z)(Zq  -  z)  Zq  -  asine 

which  can  be  integrated  by  using  standard  integral  tables. 

Burlington  (1973)  gives  the  value  of  the  integral  in  (D.52)  as 


d0 

/ - 

-  asinG 


(D.53) 


Making  the  appropriate  substitutions  of 


-1 


-a  +  ZpSinQ  +\ja^  -  Zq^cosg") 


y/ 


a^  -  z  2 


Log 


z„  -  sinG 
0 


J 


si  no  =  z/a 
cosG  =  il(z)/a 

*■  V  =  ”(^o^ 


gives  the  result  of  (D.48)  as 
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dz 

^  I(z)(Zq  -  z) 

(D.54) 

The  development  of  equation  (D.54)  was  restricted  to  values  of 
|z|  <  a.  However,  for  values  of  |z|  >  a  the  substitution  z  =  aCSC  0 
can  be  made  and  a  similar  process  repeated.  The  results  of  this  integra¬ 
tion  are  identical  to  equation  (D.54).  Consequently,  equation  (D.54) 
is  valid  for  all  values  of  z  (the  same  can  be  shown  to  be  true  for 
all  values  of  z^). 

After  much  labor  and  simplification  the  integral  of  B(z,Zq) 
was  found  to  be 


W 


Log 


-a*  +  zZq  -  I(Zq)I(z)* 
a(Zo  -  z) 


Bl(z,z-)  =  -  <  Log(z  +  z.)  -  Log(z  -  z)  +  n  \Log(zQ  -  z)  -  Log(z +  z^ 


-  "o> 


I(z) 

[z  -  z. - XI(z,z  )  -  nXI(z 

”  I(Z„)  ' 


’^o)| 

(D,55) 


where 


(zZq  -  a^  -  I(z)I(Zq))(Zq  +  a)‘ 

XI(z,z^)  =  Log 

'  (zzp  +  a^  -  I(z)I(zQ))(Zg  -  z) 
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Therefore,  integration  of  equation  (D.43)  can  be  expressed  as 

(|)(z,ZQ)sa)(z,ZQ)  =  /$(z,ZQ)dz  =  /Q(z,Zjj)dz  =  SBI(z,Zq)  +  ^BI(z,Zq) 

(D.56) 

By  integrating  equation  (D.45),  \|j(z)  can  be  expressed  in  terms  of 

equation  (D,56)  as 

ij;{z)  =  ^(z)  -  z$(z)  (D.57) 

Therefore,  in  summary,  the  functions  required  to  evaluate 
equations  (D.l)  through  (D.4)  for  the  cracked  metal  sheet  which  has  a 
stress  applied  to  the  crack  surface  (equal  in  magnitude  but  of 
opposite  sign  to  the  stress  along  the  crack  line  for  a  solid  sheet) 
are  given  by  the  equations  (D.43),  (D.45),  (D.56),  and  (D.57). 

Green's  Functions  for  Stress 

As  mentioned  previously,  the  solution  for  the  cracked  sheet  is 
obtained  from  superposition  of  the  stresses  from  the  two  loading 
conditions  shown  on  figure  D.l a  and  D.lb.  The  stress  functions 
for  these  two  loading  conditions,  which  were  developed  in  the  previous 
sections,  were  used  to  obtain  the  stresses  in  each  loading  condition. 

The  stresses  for  the  two  conditions  were  then  added  to  form  the  Green's 
function  for  stresses  for  the  loading  condition  shown  in  figure  D.lc. 
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The  stress  state  for  four  point  loads  acting  on  a  solid  sheet 
shown  on  figure  D.la  was  found  by  substituting  equation  (D.6), 

(D.7),  and  (D.8)  into  equations  (D.l)»  (D.2),  and  (0,3).  The  result 
is,  in  terms  of  the  coefficients  of  the  X  and  Y  point  loads: 

=  ReaHGl  -  G2  -  G31X  +  ReaHGlA  -  i(G2  -  G3)}Y 

Oy  =  ReaUGl  +  G2  +  G3}X  +  ReaHGlA  +  i(G2  -  G3)}Y  (D.58) 

0  =  Img  {G1  +  G2  +  G3}X  +  Img  (GIA  +  i(G2  -  G3)}Y 

xy 

where 


Real 

27r(l  +  r))tfp 


(D.59) 


(D.60) 


G2 


1 

2.11(1  +  n)t^ 


(D.61) 
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27r(l  +  n)t^  I 


(z  +  z^)^  (z-Io)^ 


(D.62) 


The  stress  state  for  the  stress  applied  to  the  crack  surface, 
as  shown  on  figure  D.lb,  was  found  by  substituting  equations  (D.43), 
(D.45),  and  (D.46)  into  equations  (D.l),  (D.2),  and  (D.3).  After 
several  algebraic  manipulations,  the  result  was  found  to  be 


=  Real{2G5  -  G7}  X  +  Real{2G6  -  G8}Y 


Oy  =  Real{2G5  +  G7}  X  +  Real{2G6  +  G8}Y 


(D.63) 


where 


=  Img  {G7)  X  +  Itng  (G8)  Y 


2Tr(l 

+  n)t 

m 

i 

27r(l 

* 

.  (z  - 

■  Z) 

2n(l 

*  ">*111 

id 

-  z) 

2tt(1 

*  ">*r.i 

B(z,Zjj)  +  B(z,Zq) 


B(z,Zq)  -  B(z,Zq) 


B'{z,Zq)  +  B'(z,Zq) 


B'(z,Zq)  -  B'(z,Zq) 


(D.64) 


(D.65) 


(D.66) 


(D.67) 
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Superimposing  the  stresses  from  the  two  loading  conditions  by 
adding  equations  (D.58)  and  (0.63)  and  taking  the  coefficients  of  the 
X  and  Y  loads  for  the  different  stresses  yielded  six  Green's 
functions  for  stress  as 


GS  =  ReaHGl  -  G2  -  G3  -  2G5  +  G7} 

XX 

(D.68) 

GS  =  ReaHGlA  -  i{G2  -  G3)  -  2G6  +  r78} 
xy 

(D.69) 

GS„^  =  ReaHGl  +  G2  +  G3  -  2G5  -  G7} 
yx 

(D.70) 

GSyy  =  ReaHGlA  +  i(G2  -  G3)  -  2G6  -  68} 

(D.71) 

(xy)x  "  +  G2  +  G3  -  G7)} 

(D.72) 

(xy)y  ~  ^  -  G3  -  G8)} 

(D.73) 

where  the  first  index  on  GS  indicates  the  stress  and  the  second 
indicates  the  load  responsible  for  it.  For  example,  indicates 

the  Oy  stress  at  point  z  due  to  a  unit  load  applied  in  the  x- 
di recti on  at  point  z^. 

The  Green's  functions  given  by  equations  (D.68)  through  (D.73) 
were  verified  with  the  use  of  a  finite  element  program  developed 
by  Y.  K,  Cheung  and  I.  P.  King  and  documented  in  Zienkiewicz's  book 
(1971).  The  finite  element  model  used  for  the  test  case  is  shown  in 
figure  D.3.  Because  of  symmetry  only  the  first  quadrant  of  the  cracked 


Fig.  D.3.  Finite  element  mesh  used  to  check  Green's  functions 
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sheet  was  modeled;  the  crack  was  simulated  by  freeing  the  nodes  in  the 
y-direction  along  the  x-axis  from  the  origin  to  x  =  1.5  inches. 

Point  loads  of  X  =  1  and  Y  =  2  were  applied  to  the  model  at  the 
point  Zq  =  2.5  +  2.5i. 

The  finite  element  results  were  compared  to  the  Green’s  functions 
results.  For  the  comparison,  the  stresses  were  calculated  along  the 
line  z  =  X  +  1.25i  by  equations  (D.68)  through  (D.73)  and  with  finite 
elements.  Finite  element  values  of  stresses  along  this  line  were 
taken  as  the  average  of  two  elements  midway  between  y-coordi nates  of 
the  nodes.  On  figure  (0.4)  the  dotted  line,  dashed  line,  and  solid 
line  indicate  a  ,  a„,  and  stresses  respectively,  obtained 

with  the  Green’s  functions  while  the  symijols  represent  stresses  obtained 
from  the  finite  element  solution.  The  comparison  was  good  and  verified 
the  Green's  functions  within  the  accuracy  of  the  finite  element  model. 

Green’s  Functions  for  Displacements 

The  displacement  field  for  the  solid  sheet  under  four  point 
loads  as  shown  on  figure  D.la  was  found  by  substituting  the  stress 
functions  $(z),  ,p(z),  and  ,|;(t)  shown  in  equation  (D.8),  (D.9), 

and  (D.IO)  respectively  into  equation  (D.4).  The  result  is 

2G(u  +  iv)  =  SFp(z,Zq)  +  SFj,(z,Iq)  (D.74) 


where 


X,  inches 


Verification  of  Green's  functions  for  stresses  in  a 
cracked  isotropic  sheet 
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(D.75) 


The  displacement  field  for  stress  on  the  crack  as  shown  in  figure 
D.lb  was  found  by  substituting  the  equations  for  the  stress  functions 
ip(t),  and  ${z)  shown  in  equations  (D.56),  (D.57),  and 
(0.43)  into  equation  (D.4).  The  result  is 


2G{u  +  1v)  =  SFq(z,Zq)  +  SFqCz.Zq) 


(D.76) 


where 


Fq(z,2q)  =  -n  BI(z,Zq)  +  BI(z,Zq)  -  (z  -  z)B(z,Zq)  (D.77) 


Superimposing  the  displacement  equations  for  the  two  preceding 
equations  yields  the  displacement  field  for  the  load  condition  shown 
on  figure  D.lc  as 


2G{u  +  iv)  =  S'Fp(z,Zq)  - 


(D.78) 


.169 
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The  Green's  functions  for  displacements  were  obtained  from 
equation  (D.79)  by  forming  coefficients  of  the  X  and  Y  loads 
for  the  u  and  v  displacements.  The  result  was 


"  S  ReaHg(z.Zo)  * 

DGxy  =  Cq  Real(i(g(z,ZQ)  -  g(z,Zg)))  (D.81) 
DGyx  '  ^0  (9(z,Zq)  +  g(z,ZQ))  (D-82) 
DGyy  =  Cq  Img  (i(g(z,ZQ)  -  gCz.z^)))  (D.83) 


4Gt  ir(l  +  n) 
m 


with 
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where  the  first  index  of  GD  indicates  the  displacement  and  the 
second  indicates  the  load  responsible  for  it. 

The  Green's  functions  for  displacements  given  in  equations 
(D.80)  through  (D.83)  were  verified  with  the  finite  element  model 
discussed  in  the  previous  section  and  shown  on  figure  D.3.  A 
comparison  of  the  displacements  calculated  from  the  Green's  functions 
and  the  finite  element  solution  is  shown  on  figure  D.5.  The  comparison 
is  made  along  the  line  z  =  x  +  2i  where  the  finite  element 
displacements  are  taken  from  the  nodal  points  of  the  model.  In  the 
figure  the  solid  and  dotted  lines  indicate  the  u  and  v  displacements 
calculated  with  the  Green's  functions  while  the  symbols  represent 
values  obtained  from  the  finite  element  solution.  The  agreement  is 
good  and  verifies  the  Green's  functions  within  the  accuracy  of  the 
finite  element  model. 


displacement,  in 


APPENDIX  E 


GREEN'S  FUNCTION  FOR  AN  UNCRACKED 
ORTHOTROPIC  SHEET 

Lekhnitskii  (1968)  gives  the  stress  and  displacements  in  an 
orthotropic  sheet  in  plane  stress  in  terms  of  two  stress  functions, 
$j(Zj)  and  $2(^2)’ 


2Real  |si  +  S2^$2'(z2)| 
2  Real|'$j'(zi)  + 

-2Real|si$i  ■  (zj  +  S2$2'{z2)| 
2ReaL  Pj^jCZi)  +  ' 
2Real|^gj$j(z^)  +  q2$2(^2^j 


)  +  P2f2’(z2)f  -  w^y  + 


)  +  q2^>2(22)>  +  w^x  +  Vq 


where  Sj  and  s^  are  the  roots  of  the  equation. 


3.  .  A  ...La  .  =  0 


(E.1) 


(E.2) 


(E.3) 


{E.4) 


(E.5) 


(E.6) 
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where  E  ,  E  ,  G  ,  and  v  are  the  elastic  moduli  and  Poisson's 
X  y  xy  xy 

ratio.  Leknitskii  (1968)  proved  that  the  roots,  s^  and  s^, 

could  not  be  purely  real  for  real  materials  but  are  purely  imaginary 

or  complex.  The  complex  roots  occur  in  conjugate  pairs  and  s^  and 

s  are  distinct  roots,  i.e.  not  complex  conjugates.  The  following 
2 


were  defined 


''xy^x 


ij  =  X  +  Sjy  Z2  =  X  +  S2y 


Pi '  r  ^ '  r  ■  '’xv’ 


(E.7) 


q  = -  (1  -  s  q 

Hi  .  p  yx  1  '  ^ 

1  y 


= -  (1  -  s 

^  s  E 


(E.8) 


Note  that  w^,  u^,  and  in  equation  (E.4)  and  (E.5)  are  rigid 

body  rotations  and  translations  respectively  and 


d<^j(Zj) 


,  $2'(z  )  = 


The  two  required  stress  functions  for  a  point  load  acting  on  a 
solid,  orthotropic  sheet  were  given  by  Lekhnitskii  (1968)  as 


$i(z,)  =  Aj,Logz^ 


(E.9) 
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where  A  and  B  satisfy  the  equations  (for  material  axis  concurrent 
c  c 

with  the  axis  of  orthotropicity) 


With  the  use  of  the  preceding  equations,  stress  functions  for 
point  loads  acting  on  a  unidirectional  boron/epoxy  composite  were 
developed.  For  this  material  the  material  constants  are 


=  0.27  X  10^  psi  E  =  3.0  x  10^  psi 

X  j 

=  0.7  X  10^  psi  =  0.019 


With  these  material  constants,  the  roots  of  equation  (E.6)  were  found 
to  be  purely  imaginary  as 
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s  =  ±0.1531  s  =  ±1.951 

1,3  2,4 

For  purely  Imaginary  roots  equations,  and  were  found  from 
equations  (E.IO)  as 

=  Cll  X  +  1  C12  Y  (E.ll) 

B  =  C21  X  +  1  C22  Y  (£.12) 

c 


C21  =  - 


S  + 

2  xy 


C22  =  - 


4irt  (s  2  -  s 
c  2  1 


Using  equations  (E.ll)  and  (E.12)  and  translating  the  origin 
so  that  the  singularity  occurs  at  the  point  z^,  equations  (E.9) 
became 


$j(Zi,Wi)  =  (Cll  X+  i  C12  Y)Log(z^  -  w^) 
$^(z^,w^)  =  (C21  X+  i  C22  Y)Log(z^  -  w^) 


(E.13) 
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where  w  =  +  s  w  =  x.  +  s  y. 

jOjO  2^2^ 

Equations  (E.13)  were  used  to  construct  a  solution  for  four 
point  loads  acting  on  a  solid  orthotropic  sheet  as  shown  in  figure  E.l. 
The  result  was 

$  (z  ,w  )  =  Cll  G8(z  ,w  )X  +  i  C12  G9(z  ,w  )Y 
111  11  11 

(E.14) 

$  (z  ,w  )  =C21G8(z  ,w  )X  +  i  C22  G9(z  ,w  )Y 
2  2  2  2  2  2  2 

where 

G8(z,w)  =  Log(z  -  w)  -  Log{z  +  w)  -  Log(z  +  w)  +  Log(z  -  w) 

(E.15) 

G9(z,w)  =  Log(z  -  w)  +  Log{z  +  w)  -  Log(z  +  w)  -  Log(z  -  w) 

(E.16) 

The  derivatives  of  equations  (E.14)  were  found  to  be 


$  '(z  ,w  )  =  Cll  G8'(z  ,w  )X  +  i  C12  G9'(z  ,w  )Y 
111  11  22 

(E.17) 

(D  '(z  ,w  )  =  C21  G8'(z  ,w  )X  +  i  C22  G9'{z  ,w  )Y 
2  2  2  1  1  2  2 


where 


nt  loads  on  orthotropic  sheet 
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G8'(z,w) 


2w  2w 

- -  +  - 


(E.18) 


G9'(z,w) 


2w 


(E.19) 


Equations  (E.13)  and  (E.17)  can  be  used  in  equations  (E.l) 
through  (E.5)  to  describe  the  stresses  and  displacements  in  the 
orthotropic  sheet. 

Green's  Functions  for  Stresses 

The  Green's  functions  for  stresses  were  developed  by 
substituting  equations  (E.17)  into  equations  (E.l),  (E.2),  and  (E.3) 
to  determine  the  stress  state  as 


“x  =  ''SxxX  * 


(E.20) 


“xy  - 


where  the  Green's  functions  are  given  by 
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2Realls^^  Cll  G8'  (z^,  w^)  +  C21  G8' 


.«,)j 


■  ■' 

\<  i<  s  ‘ 
1 

.(.  (. 


2Real|  ijs  ^  C12  G9'  (z^.w^)  +  C22G9'  (z^.w^ll  (E.24) 


=  2Real|  Cll  G8'  (z^,w^)  +  C21  G8'  (z^.w^)|  (E.25) 


2Real|  i  j  C12  G9'  (z^,w^)  +  C22  G9'  (z 


*(xy)x 


’(xy)y 


■2Realls  Cll  G8'  (z  ,w  )  +  s  C21  G8'  (z 
J  1  112 


Z  ,W  )] 
2  2  J 


(E.26) 


(E.27) 


■2Real<  i<s  C12  G9'  (z  ,w  )  +  s  C22  G9'  (z 
111  112 


-I 


(E.28) 


To  verify  these  Green's  functions,  stresses  computed  with  equa¬ 
tions  (D.20)  through  (D.22)  were  compared  to  finite  element  results. 
The  model  used  for  the  finite  element  solution  is  identical  to  that 
shown  in  Appendix  D  on  figure  D.3  except  that  no  nodes  were  freed 
along  the  x-axis  to  simulate  a  crack.  The  finite  element  program  used 
is  documented  by  Zienkiewicz  (1971)  as  mentioned  in  the  previous 
appendix.  However,  the  program  as  given  by  Zienkiewicz  does  not  have 
orthotropic  capability.  Therefore,  it  was  modified  by  introducing 
the  orthotropic  stiffness  matrix  for  plane  stress  given  by 
Zienkiewicz  (1971)  as 
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[D]  = 


1  -  nv 


yx 


nv 


yx 


where  n 


and  m  = 


yx 


nv  0 

yx 

1  0 

0  M(l  -  nVyJ^^) 

in  place  of  the  isotropic 


stiffness  matrix  in  the  program. 

The  comparisons  between  the  stresses  along  the  line  z  =  x  + 

1.25i  are  shown  on  figure  E.2  on  which  the  dotted  line,  dashed  line, 

and  solid  line  indicate  the  a  ,  a  ,  and  a  stresses  respectively 

X  y  xy 

from  the  Green's  functions.  The  symbols  on  the  figure  indicate  the 
same  stresses  obtained  from  the  finite  element  solution.  The 
comparison  was  good  and  verified,  within  the  accuracy  of  the  finite 
element  model,  the  accuracy  of  the  Green's  functions  for  stresses. 


Green's  Functions  for  Displacements 


The  Green's  functions  for  displacements  were  developed  by 
substituting  equations  (E.14)  into  equations  (E.4)  and  (E.5).  Because 
of  the  symmetry  of  the  loads  the  rigid  body  rotation,  w^,  and  the 
rigid  body  translations,  u^  and  v^,  are  zero  in  equations  (E.4)  and 
(E.5).  The  results  are 


(E.29) 


V  •  HD^^X  . 


(E.30) 


stress,  psi 


Fig.  E.2.  Verification  of  Green's  functions  for  stresses  in 
an  uncracked  orthotropic  sheet 
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where  the  Green's  functions  are  given  by 


=  2Real  {p  Cll  G8(z  ,w  )  +  p  C21  G8(z  ,w  ) 

XX  1  1  112  22 


=  2Real 


C12  G9(z  ,w  )  +  p  C22  G9(z 
112 


HD,„  =  2Real<  q  Cll  G8(z  ,w  )  +  q  C21  G8(z  ,w  ) 


HD. 


yy 


=  2Real  i  |q^ 


1  1  2 


2  2 


C21  G9{z  ,w  )  +  q  C22  G9{z  ,w 
112  2  2 


'] 

-’ll 

'] 


(E.31) 

(E.32) 

(E.33) 

(E.34) 


With  the  use  of  the  finite  element  model  shown  in  figure  D.3, 
equations  (E.31)  through  {E.34)  were  verified  by  comparing  the  displace¬ 
ments  along  the  line  z  =  x  +  2i  calculated  with  equation  (E.29) 
and  (E.30)  shown  on  dotted  and  dashed  line,  respectively,  on  figure 
E.3  with  the  displacements  calculated  with  the  finite  element  program. 
The  comparison  is  good  and  within  the  accuracy  of  the  finite  element 
model  verifies  the  Green's  functions  for  displacements  in  the 
uncracked  orthotropic  sheet. 


displacements,  in 


APPENDIX  F 


COMPUTER  PROGRAM  TO  PREDICT  CRACK 
AND  DEBOND  GROWTH 

THg  analysis  devGloped  in  ChaptGrs  IV  and  V  to  predict  crack 
growth  in  reinforcGd  systems  requires  the  use  of  a  digital  computer 
to  perform  numerical  integration  and  solve  large  systems  of  simultaneous 
equations.  The  analysis  was  programmed  in  FORTRAN  IV  for  use  on  the 
NASA  Langley  Research  Center  (LRC)  CDC  6600  computing  system.  However, 
no  special  system  routines  were  used  in  the  program,  and  the  program 
should  be  usable  on  any  computing  system  that  uses  a  FORTRAN  IV 
compiler.  With  the  exception  of  the  Gaussian  elimination  subroutine 
SIMQ,  the  program  is  all  original  code. 

For  economical  reasons  the  user  should  be  aware  of  the  central 
processing  time  (CP)  and  central  memory  (core)  required  to  execute 
the  program.  The  CP  time  is  primarily  a  function  of  the  numerical 
integration  of  the  Green's  functions.  Several  integrations  are 
required  for  each  shear  element,  i.e.  elements  of  region  B  shown  on 
figure  8.  Consequently,  the  CP  time  requirement  is  a  function  of  the 
number  of  shear  elements  and  can  vary  from  a  few  to  several  thousand 
seconds  depending  upon  the  number  of  elements.  The  core  requirements 
are  also  a  function  of  the  number  of  shear  elements.  Each  element 
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contains  two  unknown  shear  stresses,  and  The  number  of 

rows  and  columns  of  the  square,  fully  populated  coefficient  matrix 
used  to  solve  for  these  stresses  is  two  times  the  number  of  elements. 
In  the  program  the  dimensions  are  set  for  fifty  shear  elements 
(100  unknowns),  but  the  program  can  accomodate  more  elements  if  the 
array  Z9(10,000)  in  the  program  is  enlarged.  For  fifty  elements  the 
core  requirement  is  137K  octal. 

Data  are  read  in  the  program  via  a  NAMELIST  with  the  following 
definitions: 

E3  modulus  of  the  cracked  sheet 
T2  thickness  of  the  cracked  sheet 
VI  Poisson's  ratio  of  the  cracked  sheet 
E4  modulus  of  the  composite  sheet  in  the  y-direction 
(loading-axis) 

E5  modulus  of  the  composite  sheet  in  the  x-direction 
V4  Poisson's  ratio  of  the  composite  sheet  for  a  load  applied 
in  the  x-direction 

GC  shear  modulus  of  the  cracked  sheet 
T4  thickness  of  the  reinforcement  sheet 

TAD  thickness  of  the  adhesive  layer 

GAD  initial  effective  shear  modulus  of  the  adhesive 
GAD2  secondary  effective  shear  modulus  of  the  adhesive 

SYIELD  yield  stress  of  the  bulk  adhesive  in  uniaxial  tension 

A1  initial  crack  length  in  metal  sheet 

F  initial  aspect  ratio  of  the  debond  ellipse 
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NC  number  of  columns  of  interlaminar  sheat  elements 

NR  number  of  rows  of  interlaminar  shear  elements 

FCYC  final  number  of  applied  load  cycles 
S  remote  applied  stress  applied  to  the  reinforced  system 
in  the  y-di recti on 

As  an  example  a  sample  run  was  made  for  one  iteration.  Figure 
F.la  shows  the  model  prior  to  the  iteration  and  figure  F.lb  shows 
it  after  the  iteration.  The  sample  NAMELIST  input  for  the  run  was 


as  follows: 


25.500  psi  5  ,  25.500  psi 
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o 


(a)  Model  before  iteration  (b)  Model  after  iteration  (12  cycles) 
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SAMPLE  INPUT 

$INPU  E3=l . OE+07 .T2= . 1 56 , VI= . 30, E4=3 . OE+07 , E5=3 . OE+06 , V4= . 02 ,GC=7 . OE+05 , 
TAD=.004,GAD=65000. ,GAD2=36000. ,SYIELD=4200. ,A1=1 .00,F=.001  ,NC=5,NR=3, 
FCYC=10000. ,S=25500. ,T4=.0208, 

$ 


SAMPLE  OUTPUT 


MODULUS  OF  THE  CRACKED  SHEET  .107E+08 
THICKNESS  OF  THE  CRACKED  SHEET  .156E+00 
POISSONS  RATIO  FOR  THE  CRACKED  SHEET  .300E+00 
COMP  MODULUS  PARALLEL  TO  LOAD  .300E+08 
COMP  MODULUS  TRANSV  TO  LOAD  .300E+07 
SHEAR  MODULUS  OF  COMPOSITE  .700E+06 
COMP  POISSONS  RATIO  TRANSV  TO  LOAD  AXIS  .200E-01 
THICKNESS  OF  THE  REINFORCEMENT  SHEET  .280E-01 
MINOR  AXIS  HEIGHT  IN  PER  CENT  CRACK  LENGTH  .lOOE-02 
REMOTE  APPLIED  STRESS  .255E+05 
NUMBER  OF  COLUMNS  FOR  BOUNDARY  POINTS  5 
NUMBER  OF  ROWS  FOR  BOUNDARY  POINTS  3 
FINAL  NUMBER  OF  APPLIED  LOAD  CYCLES  .lOOE+05 
INITIAL  CRACK  LENGTH  .lOOE+01 
ADHESIVE  THICKNESS  .400E-02 
SHEAR  MODULUS  OF  THE  ADHESIVE  .650E+05 
SECONDARY  SHEAR  MODULUS  AFTER  YIELDING  .360E+05 
YIELD  STRESS  OF  MODULUS  IN  UNIAXIAL  TENSION  .420E+04 


METAL  SIGM-X=  .269E+03  SIGM-Y=  .200E+05 

COMPOSITE  SIGM-X=  -.150E+04  SIGM-Y=  .559E+05 

X  AND  Y  DIMENSIONS  OF  INTEGRATION  AREA  ARE  .990E+00  .746E+00 
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NODE  NUMBER 


X-COR 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


.990E-01 

.297E+00 

.495E+00 

.693E+00 

.891E+00 

.990E-01 

.297E+00 

.495E+00 

.693E+00 

.891E+00 

.990E-01 

.297E+00 

.495E+00 

.693E+00 

.891E+00 


14  ELEMENTS  HAVE  YIELDED  AT  .189E+05 
YIELDED  ELEMENTS 

3  2  14 

8  7  10  6 


THE  YIELD  MACROSCOPIC  STRESS  IS  .716E+04 


NODE  DEBOND  HEIGHT 


1 

.125 

2 

.121 

3 

.113 

4 

.099 

5 

.076 

6 

.373 

7 

.365 

8 

.349 

9 

.323 

10 

.285 

11 

.622 

12 

.613 

13 

.593 

14 

.563 

15 

.521 

ELASTIC 

X-COEFFICIENT  Y 

.300E+01 

.279E+02 

.lOOE+03 

.327E+03 

.lllE+04 

.346E+02 

.130E+03 

.328E+03 

.806E+03 

.191E+04 

.916E+02 

,293E+03 

.560E+03 

.983E+03 

.180E+04 


Y-COR 

.125E+00 
.121E+00 
.113E+00 
.988E-01 
.765E-01 
.373E+00 
.365E+00 
.349E+00 
.323E+00 
. 285E+00 
.622E+00 
.613E+00 
.593E+00 
.563E+00 
.521E+00 


-COEFFICIENT 

.849E+04 

.859E+04 

.864E+04 

.847E+04 

.667E+04 

.256E+04 

.266E+04 

.276E+04 

.276E+04 

.220E+04 

.117E+04 

.119E+04 

.117E+04 

.108E+04 

.807E+03 
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WITH  PLASTICITY 


NODE 

DEBOND  HEIGHT 

X-COEFFICIENT 

Y-COEFFICIENT 

1 

.125 

.207E+01 

.708E+04 

2 

.212 

.213E+02 

.709E+04 

3 

.113 

.795E+02 

.704E+04 

4 

.099 

.245E+03 

.679E+04 

5 

.076 

.805E+03 

.536E+04 

6 

.373 

.438E+02 

.312E+04 

7 

.365 

. 1 56E+03 

.321E+04 

8 

.349 

.371E+03 

.325E+04 

9 

.323 

.864E+03 

.315E+04 

10 

.285 

.191E+04 

.244E+04 

11 

.622 

.105E+03 

.165E+04 

12 

.613 

.333E+03 

.165E+04 

13 

.593 

.635E+03 

.161E+04 

14 

.563 

.110E+04 

.146E+04 

15 

.521 

.189E+04 

.112E+04 

K-BOND 

K-UNSTIFFEND 

K-STIFFENED 

K- FACTOR 

-.108E+05 

.200E+05 

.919E+04 

.459E+00 

STRESS  AND  STRAIN 

IN  THE  METAL 

NODE 

SIG-1 

SIG-2 

SIG-12 

1 

.681E+04 

.580E+04 

-.472E+02 

6 

. 340E+04 

.121E+05 

-,103E+03 

11 

.137E+04 

.154E+05 

-.146E+03 

NODE 

EPS-1 

EPS-2 

EPS-12 

1 

-.799E-03 

.733E-03 

-.115E-04 

6 

-.657E-03 

.123E-02 

-.250E-04 

11 

.559E-03 

.148E-02 

-.354E-04 
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STRESS  AND  STRAIN  IN  THE  COMPOSITE 


NODE 

SIG-1 

SIG-2 

SIG-12 

1 

.838E+04 

.105E+06 

.742E+03 

6 

.897E+04 

.728E  05 

.252E+03 

11 

.758E+04 

.568E+05 

-.301E+03 

NODE 

EPS-1 

EPS-2 

EPS-12 

1 

,272E-02 

.350E-02 

.674E-04 

6 

.294E-02 

.242E-02 

.360E-04 

11 

.249E-02 

.189E-02 

-.430E-04 

ENERGY  RELEASE  ON 

MINOR  AXIS 

IS 

.597E+01 

APPLIED  CYCLES= 

. 990E+01 

DA/DN  BEFORE  CYCLE 

INCREMENT 

IS 

.181E-04 

DB/DN  BEFORE  CYCLE 

INCREMENT 

IS 

.202E-01 

INCREMENT  CYCL  =  .990E+01  CRACK, LENGTH=  .lOOE+01 


MAX  DEBOND  HEIGHT=  .201E+00 


ooooooooooooooooooo 
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IROGRAM  COLLS(  INPUT , OUTPUT, TAPE4) 

C 

C  THIS  PROGRAM  CALCULATES  CRACK  AND  DEBONL  GROWTH  IN  A 
C  REINPORCED  SYSTEM  COMPOSED  OF  A  CRACKED  METAL 
C  SHEET  REINFORCED  V/ITH  A  UNIDIRECTIONAL  BORON/EPOXY 
C  SHEET. 

C 

COMMON /R0T/Z9(10000),D(  100  ),NALF 

00MM0N/B0ND2/NE  ,NL  ,NT  ,XD(  100  )  ,YC(  100  )  ,XA(  100  )  ,YA(  100  ) , 
INOP(IOO) 

DIMENSION  FTNC(50),EXNCY(50),FA1(50),FB1(50  ) 

COMMON /ADHES/TAD ,  GAD 

COMMON  /TO  P/E3 ,  T2 ,  V 1 ,  SMX  ,  S^IY ,  G ,  CONS ,  Q ,  A1 ,  SCON 
DIMENSION  DR(  100  )  ,GPR(  100  )  ,IiPR(  100  ) 
COMMai/XLIMIT/DC(2,51  ),DB(6,51  ),EP 
DIMENSION  CM(  3 , 3 )  ,SIG ( 3) ,  STRAIN(  3  ) .  DD(  100  ) 

DIMENSION  SM(3,3),CI-IM(3,3),CMC(3,3),SrRESSM(3),SrRESSC 

DIMENSION  TSIGM(3),TSIGMT(3),ayiM(3,5),SMC(  3,3) 
DIMENSION  TOTG(  50  ) ,  SIGMT(  3) ,  TSTRAIN(  3 ) ,  STRESS (2 , 3 , 50 ) 
1,STRANN(2,3,50) 

COMMON /BO T/E4,  T4,  V4,  SCX  ,SCY  ,GS  ,CONB  .Q1,  GC  ,E5 
COMMON/CTOL/TOL,NC  ,NR  ,TX  , TY  ,NBC ( 100  )  ,LBC 
COMM®  /BO ND/F ,  P 1 ,  P2 ,  XKF  , XKUN S , XKS T  IF 
COMM® /lEC  AY/EE  ,ALPHA1  ,ALPHA2 
coMPiEX  z,ci,xk: 

EXTERNAL  XK 
CI=CMPLX(0.  ,1.0) 

NAMELIST/INPU/E3,  T2, F , S, TEST  ,T4,  E4,  V4,  V1 , E5,  GC  ,  TAD, GAD 
1,NR,NC  ,FGYG, 

1A1,GAD2,SYIEED 

E3-M0DULUS  OF  THE  CRACKED  SHEET 
T 2- THICKNESS  OF  THE  CRACKED  SHEET 
VI- IS  THE  POISSONS  RATIO  FOR  THE  CRACKED  SHEET 
GC -SHEAR  MODULUS  OF  THE  REIl'IFORCMENT  SHEET 
E4-MDDULUS  OF  THE  REINFORCEMENT  SHEET  IN  THE  LOADING 
DIRECTION 

E3-M0DULUS  OF  THE  REINFORCEMENT  SHEET  TRANSVERSE  TO 
THE  LOADING  AXIS 

V4  IS  THE  POISS®S  RATIO  FOR  THE  BOTTOM  SHEET  FOR  A 
TRANSVERSE  LOAD 

T4-aHICKNESS  OF  THE  REINFO RCEI-IENT  SHEET 
F-MINOR  AXIS  HEIGHT  IN  PERCENT  CRACK  LENGTH 
S  -REMOTE  APPLIED  STRESS 
NH-NUMBER  OF  ROV/S  OF  BOUNDARY  POINTS 
NC-NUMBER  OP  COLUMNS  OF  BOUNDARY  POINTS 
FCYC-HNAL  NUMBER  OF  LOAD  CYCLES 
AI-naTIAL  CPvACK  LMGTH 


K=0 


o  o  o 
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1000  OONTinUE 
READ  INPU 
MX=RR>^NC 
M=MX 
N=MX 
NALF=M 

OUTPUT  TOE  IfIPUT 

PRINT  1,E5 
PRINT  2,T2 
HUNT  20, VI 
HINT  18,  E4 
HINT  26, E5 
PRINT  27, GO 
HINT  21,  V4 
PRINT  19,  T  4 
HINT  3,F 
PRINT  7,S 
PRINT  30, HC 
HINT  31,  NR 
HINT  32,FCYC 
HINT  34, A1 
HINT  24, TAD 
HINT  25,  gad 
HINT  5,GAD2 
HINT  e.SYIEED 
P1=2. 

P2=2. 

T0L=.001 

FF=F 

N=2^N 

M=2*M 

NSQ=N*N 

IF(N.GT.100  )  500, 501 

500  HINT  502, M,N 

502  FCRMAT(  //*  ERRORS  —  M  OR  N  EXCEEDS  DIMENSIONS  BOUNDS 
1M=*I10*N=*I10//) 

GO  TO  1001 

501  CONTINUE 
Q=(3.-V1)/(1.+V1) 
q1=(3.-V4)/(1.+V4) 

G=E3/U.'*^h  .+V1)) 

C0NS=1./(12.566*T2  *(l-tQ)*G) 

SC0N=1  ./(6.2852*(  1+Q  )*T2) 

GB=E4/(2.*(  1.+V4)) 

C0NB=1./(12.56  6*T4  *(14Q1)*GB) 

V5=V4*E5/E4 

ES=2. 7182818 

ALPHA.  1=GAD*(  1 .  /  (T  2*E3  >+  1 .  /  (T  4*E4 ))  /TA.D 
ALPHA  1=SQRT  ( AU^Ek  1 ) 
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ALPHA2=GAD2^(  1 .  / (T  2*E  3 )+  1 .  / (T  4*E 4 ))  /^D 
ALPHA  2=SQRT(  ALPHA  2 ) 

TY=-ALOG(  .05  )/ALPHA2 
C 

C  CALCULATE  THE  REMOTE  STRESSES 
C 

C  ZERO  ALL  COMPHANCE  MATRICES 
C 

DO  50  1=1,3 
DO  50  J=1,3 
CI'1(I,J)  =  0. 

CMM(I,J)  =  0. 

cr4c(i,j)=o. 

50  CON T  HUE 

c 

C  CALCULATE  COr-ULIAlICE  EOR  METAL  SHEET 
C 

a'IM(l,l)=E3/(l.-V1*^-2) 
aiM(1,2)=V1^-CMM(l,  1  ) 

CMM(2,1)=CMM(1 ,2} 

CMM(2,2)  =  CI'M(1 ,1) 

CMM(3,5)=G 

C 

C  GENERATE  THE  STIFF  MATRIX  FOR  THE  METAL  SHEET 

C 

CALL  INYER  ( G-IM ,  ,  DD ) 

C 

C  CALCULATE  Ca-IPLIANCE  FOR  C®POSITE  3IEET 
C 

:{N=E5/E4 

.’(r-I=GC/E5 

CC=E4/(l.-M*V4**-2) 

CMC(1  , 1}=:0H*CC 
CMC(1,2)=:{N*V4*CC 
CHC(2,1)=CIC(1 ,2) 

C1/[C(2 1  2)  =  GG 

CI'IC  ( 3 !  3 ( 1 .  -XN*V4^  2  )*CC 

c 

C  GEIIERATE  OHE  STIFF  I4ATRIX  FOR  THE  COMPOSITE  SHEET 

C 

CALL  INVER(CMC,3f'IC  ,DD) 

C 

C  GENERATE  MACRO  STIFFNESS  mTRIX 
C 

DO  51  1=1,3 
DO  51  J=1,3 

CI>:  ( I ,  J)  =  cm  ( I ,  J )  *  T  2+CM  G  ( I ,  J )  -X  T 4 
51  CONTHUE 
C 

CALL  IN¥ER(  CM ,  SM  ,  ED ) 


oooo  ooo  ooo  ooo 
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3IG(1)=0.0 
SIG(2)=S*(T2+T4) 

SIG(3)=0.0 

CALL  MULT(SM,SIG, STRAIN) 

400  K)RMAr(3S11  .3) 

CAICULATE  STRESSES  IN  EACH  LAYER 

CALL  MULT (CMM. STRAIN, STRESSM) 

ST'IX^STRESSMd  ) 

3-iY=STRESSM(  2) 

CALL  MUIff  ( CMC ,  STRA  IN , STRE S^X! ) 

SCX=STRESSC(1) 

3CY=STRESSC(2) 

PRINT  300,SKX,3^IY 

300  F0RI4AT(/*  META.L  SIGM-X=*E  10 . 3-^  SIGM-Y=*E10.3) 

PRINT  30 1,  sex, SC Y 

301  P0RKAT(/*  COMPOSITE  SIGM-X=*E  10 . 5"'  SIGM-Y=*E  10. 3/ ) 

SET  Ur'  LOOPS  FOR  RIGHT  HAND  VECTOR  AND  COEFFICIENT  MATRIX 

T1=T2 
CALL  CPAR 

LOOP  ON  INCRSraWTS  OP  LOAD  CYCLES 

TNC=0.0 
K0UNT=0 
1004  CONTINUE 

IF(  KCUNT.NE  .0  )PRINT  551  ,K0UNT 
551  F0Rr'AT(//3OX*  XX:'XXXX:o:XX  njCREI-NNT  FOR  LOAD  CYCLES^- 

i*is->Y5  xxx:<xxxxx^7) 

niv_  oq-x-A  1 
PRINT  l6\T)C,aY 

16  P0RKAT(/*  X  ANDY  DIMINSiaTS  OF  WTEGRATION  AREA  ARE 
12E11  .3/) 

E0U!]T=K0UNT+1 
MSAVE=I' 

CALL  GRID 
CALL  F0Ri:(N,M) 

MAPS  RIGHT  HAI'JD  SIDE  A  UNIT  VECTOR  FOR  USE  V/ITH  PLASTIC 
AITAIYSIS 

DO  1401  1=1,  N 
C  SAVE  UNIT  RIGHT  HAND  VECTOR 
IE(I)=D(I)/S 
D(I).=D(I)/S 
1401  CONTINUE 
REWBID  4 

1130  F0RrNiT(l1E10.5) 
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C 

C  STORE  OOEFITCIEITT  MATRIX  ON  TAPE 
C 

’.VRITE(4, 1130)(Z9(I)  ,I=1,1®Q) 

PEV/IIID  4 

CALL  SIMQ(Z9,D,N,DiD) 

C 

C  RETRIVE  OOEPPICISNT  MATRIX 
C 

READ(4, 1 150)(Z9(I)  ,I=1,1^Q) 

C 

C  SAVE  ELASTIC  SOLUTiai 
C 

DO  1402  1=1,  IT 
DD(I)=D(I)-!<S 
1402  COITTINUE 
C 

C  GALL  TEE  H;/iSTIC  SUBROUTINE  TOEIITD  INTERL/MTAR  STRESSES 
C  AFTER  YIELDING 

C 

CALL  rIA3riC(SYIELD,GAD2,S,IR) 

NALE=K/2 
PRINT  17 

17  FORMATC  /70  X’! NI TH  PLAS T ICITY*29 TP^ELASTIC*  ) 

PRINT  11 
DO  104  I=1,mLP 

it:=i+iialp 

X=DC(1,I) 

Y=DC(2,I) 

PRINT  12,X,I,Y,D(I)  ,D(IK),IE)(I)  ,nD(l+NALP) 

104  CONTBTUE 
C 

:ffiT0T=0. 

DO  105  J=1,1TALP 
:<=DC(i,j) 

Y=DC(2,J) 

Z=X+CI^Y 

CALL  XLTTG(Z,J,A1,S11  ,S12,S21,S22,:«,5) 
xk;i=sii+si2 

:iIC2=-(S21+S22) 

XKTOT=xxTOT+:a:i 

105  CCaiTINUE 

c 

C  PRINT  OHE  STRESS  INTENSITIES 

c 

J<KU1T3=SITY>:  A1  .5 
.{KSTIP=Xi:UIIS+XKTOT 
/3[F=XESTIF/XvUNS 
PRINT  14 

i-RIIZT  15  ,  XICTOT ,RXUNS,XKSTI?,XKP 
100  CONTEIUE  ■ 


ooo  ooo  ooo  ooo  oooo  oooo 
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IhKSAYE 


CALCUIATE  STRAINS  AND  STRESSES  AT  COORDINATE  POINTS  IN 
THE  ADHHRaiDS 

DO  1220  I=1,I>K,11C 
X=DC(1,I) 

Y=DC(2,I) 

z=x+cr^-Y 

DO  12  26  KTYP3=1,2 

CALL  ysRi(  z ,  sa ,  ^:y  ,  s:cxy  ,  six ,  sly  ,  slxy  ,mx  ,  ktype  ) 

FIRST  THREE  TERNS  ARE  FROM  X-L6AD3  THE  LATTER  THREE  FROM 
Y-IjOADS 

TSIGM(1)=SXX+SYX 
TSia-l(2)=SXY+SYY 
TSIGM(  3)=SXmSYXY 

CALCULATE  THE  STRESSES  IN  THE  COMPOSITE  SHEET 

CALL  RESIGH(Z,S11  ,S22  ,S12  , KTYPE, SrRSSSH,STRESSC) 

SUPERIMPOSE  STRESSES 

TSI  GMT  ( 1  )=  T  SIGM  (  1  )+S  1 1 
TSIGMT(2)=TSIGM(  2}+S22 
TSIGMT(3)=TSIGK(3)+S12  ^ 

CAI;CULATS  TIDi]  CORRESPONDING  STRAINS 

IF(ETLPE.E3  .1  )  1221, 1222 

1221  CALL  I-ULT(SM!,TSIGMT,TSTRAIN) 

GO  TO  1224 

1222  CALL  r4ULT(3-IC  ,T3IGMT,TSTRAm) 

1224  CONTH’IUE 


STORE  ALL  VATUES 
DO  1225  Ii:=1,5  ^  , 

3TRES3(NTY1'E,IK  ,I)  =  TSIGMT(IK) 

STR  ANN ( NT YI'E ,  DC  ,  I)  =  T  STIU IN ( m ) 

1225  CCNTIIIUS 
12  26  OONTDIUE 
1220  CONTINUE 
PRINT  1227 

12  27  FORMATC/*  STRESS  AND  STRAIN  IN  THE  METAL  */) 

PRINT  1228 

1228  PCR!"iAT(  2X>f  N0IE*  5:C?<-3IG-1  *8>?^SIG-2*8X- SIG-12  *8X^EPS-1  *8X 
1*SPS-2*8X*SrS-12*) 

PRINT  1229,  (I,  (  STRESSd  ,K,I)  ,K=1,5) ,  (STRANNd  ,K,I)  , 
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1K=1,5),I=1,MX,NC) 

1229  F0RI'IAT(I5,  5XE10.3 ,3XE10.5,3XE10.3,4XE10.3,3XE  10.5,5XE1 

10.3) 

PRINT  1230 

1230  P0RIiAT(/*  STRESS  AND  STRAIN  IN  THE  COMPOSITE*/) 

PRINT  1228 

PRINT  1229,(1  ,(STRESS(2,K,I)  ,Ife1,3) ,  (STRANN(2,K,I)  ,E=1 

1.3) , 1=1, MX, NC) 

SUM=0. 

DX=TX/NC 

DO  1232  I=1,MX,NC 

dx=tyAir 

SUM=SUH+DY*  ( STRANN(2 ,2,1)  -STRANN(  1,2,1))  *D(MX+l) 

1232  CON T  HUE 

PRINT  1251,  SUM 

1251  K>RI'1AT(/*  ENERGY  RELEASE  ON  MINOR  AXIS  IS*E11.3/) 

CALL  XINC(  SUM  ,XNCY  ,IAS  ,IPS ) 

CPR(KOUNT)=DAS 

DPR(KOUNT)=DES 

IF(XKSTIF. (35)  .56000  .)G0  TO  1005 
TNC=TNC+XNCy 

IF(aiTC.GE.FCyC)GO  TO  1003 

B1=P*A1 

HINT  42, ®C 

42  FORMAT( /*  APPLIED  CYCLES=*E10 . 3) 

HINT  555,DAS,DPS 

555  F0RMT(*  lA/m  BEFORE  CYCLE  INCREMENT  IS*E11  .3/ 

1*  DB/M  BEFORE  CYCLE  INCREMENT  IS*E  11  .3/) 

PRINT  40,XNCY,A1,B1,P1,P2 

1*  MAXIMUM  DEBOl®  HEIGHT=*E10.3*  P1=^E10.3*  P2=*E10.3/) 
40  F0RMAT(  5X* INCREMENT  CYCL=*E10 . 3*  CRACK,LENGTH=*E  10 . 3 
PTNC(KOMT)=TNC 
FXNCY(KOltTT)=XKCY 
FA1(K0UNT)=A1 
FB1(K0UNT)=B1 
IF(A1.GT.2)G0  TO  1001 
GO  TO  1004 

1  P0m>IAT(20]P^M0DULUS  OP  THE  CRACKED  SHEET 

1  *,E12.3) 

2  FORMAT  (2O0P^THICKNESS  OF  THE  CRACKED  SHEET 

1  *,E12.3) 

3  F0H'IAT(20X?:MIN0R  AXIS  HEIGHT  IN  PERCENT  CRACK  LENGTH 

1  *,E12.3) 

4  FORMAT ( 20 X- INCREMENTS  TO  MXIUIOT  CRACK  LENGTH 
1  ->^,112) 

7  FORMAT  ( 20 REMOTE  APPLIED  STRESS 

1  *,E12.3) 

8  F0RMAT(2OX*INITrAL  CRACK  LENGTH  BEFORE  CYCLING 

1  *E12.3) 

9  FORI^AT(/25X^ INCREMENT  *I5,4X*CRACK  LENGTH  *F10.3/) 

10  FORRUiT(5E12.3) 
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11  PORI'IATCSOPLOCATION  node  debond  HEIGHT*6X 

1*X-COEFEICIENT*8}P*-y-COEFPICIENT*7X*X-COEFEICIENT* 

27)P^Y-COEFFrCIENT*/) 

12  F0RI^IA!P(22XF5.3,9XI2,10XP5.3,4(10XE10.3)) 

14  FORI'IAK  /23M-B0ND*5iP^K-UNSTIFFEND*5X>^K-STIFFENED*5X 
1*K-FACT0R*/) 

15  FOR1'IAT(20X,3(E10.3,5X),E10.3/) 

5  K)RMAT(  20 X^ SECONDARY  SIEAR  MODULUS  AFTER  YIELDING 

1  *E12.3) 

6  K)RMAT(  20 YIELD  STRESS  OF  MODULUS  IN  UNIAXIAL  TENSION 

1  *E12.3) 

18  F0RMAT(  20X^COMP  MODULUS  PARALLEL  TO  LOAD 

1  *,E12.3) 

19  F0RKAT(  20X*TniCKNESS  OF  THE  REINFOBCEMENT  SHEET 

1  *,E12.3) 

20  FORM AT(  20 POISSONS  RATIO  FOR  THE  CRACKED  SHEET 

1  *,E12.3) 

21  FORMAT(  20X*C0MP  POISSONS  RATIO  TRANSV  TO  LOAD  AXIS 

1  *,E12.3) 

22  F0RI4AT(  20X*NUMBER  OF  ELHPSES  IN  GRID 
1  *112) 

23  F0RMAT(  20X5^NUI-IEER  OF  HNES  IN  GRID 

1  *112) 

24  F0RMAT(2OX^ADHESIVE  THICKNESS 

1  *,E12.3) 

25  FORM AT(  20 X^ SHEAR  MODULUS  OF  THE  ADHESIVE 

1  *,E12.3) 

26  FORM AT(  20 X* COMP  MODULUS  TRANSV  TO  LOAD 

1  *  E 12  •  3 ) 

27  FORM AT(  20 X^ SHEAR  MODULUS  OF  COMPOSITE 

1  *,E12.3) 

30  FORMAT(  20X*NUMEER  OF  COLUI-mE  FOR  BOUNDARY  POINTS 

1  *,112) 

31  F0RMAT(  20 X* NUMBER  OF  ROV/S  FOR  BOUNDARY  POINTS 

1  *,112) 

33  FORMAT ( 20 X<-F INAL  I'lUMBER  OF  APPLIED  LOAD  CYCLES 

1  *,E12.3) 

34  F0RMAT(2OX^INITIAL  CRACK  LENGTH 

1  *,E12.3) 

1005  PRINT  58,TNC 

38  K)RMAT(5X* SPECIMEN  FAILED  BEFORE*E  10.3 *CYCLES*/) 

GO  TO  1001 

1003  PRINT  39,PCYC 

39  PORMAT(  5X*  THE  SE  ECIFIED  NUl^lBER  OF  LOAD  CYCLES  MS* 

1*  BEEN  I®T*E11  .3) 

2002  FORI«AT(/) 

1001  CONTINUE 
M=M/2 
N=N/2 

L0UNT=E0UNT-1 
PRINT  540 


oooooo  ooo  ooo 
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540  F0RHAT(/15X>^DSLTA-If«- 11  TOTAL-11^  18X'^-A-'MS]e^DA/n?!<  1SX*B^^ 

1  16X'*DB/M-^7) 

PRII1T  556,  (PXIICY(I)  ,PTUC  (l)  ,PA  1  (I )  ,aPR(  I  ),PB  1  (I)  ,IJPR(  I 
1),I=1,L0UNT) 

536  PORE AT(  6S 20.3  ) 

IP  (TEST. IIS  .0.)  GO  TO  1000 
EilD 

SUBROUTINE  IIIVER ( 3,^,1)) 


TfllS  ROUTINE  Iin^ERTS  A  5X3  EA.TRIX 

DinaiSIOi:  A(5,5),x\I(5,3),AG(  3, 3), 3(3, 3) 
I)=S(1,1)*S(2.2)^S(5,3)+S(1,2)*S(2,3)'"S(3,1)+S(1,3)* 
1S(2,1)*S(3,2) 

2-(S(  1 .3)''3(3.1)->^S(2,2)+S(1,2)*S(2,1)*S(3,3)+S(  1,1  )* 

AG(1,2)=-S(2,1)*S(3,3)+S(2,3)*S(3,1) 

AG(1,3)  =  S(2,1)^S(3,2)-S(2,2>S(3,1)^ 
AG(2,1)=-S(1,2)^S(3,3)+S(1,5)-='S(3,2) 
AG(2,2)=3(1,1)^S(5,3)-S(1,3)*S(3,1) 
AG(2,3)f“S(1,l)*S(3,2)+S(1,2)^S(3,l) 

AG(5.1  =S(1,2)-"S(2,3)-S(1,3)^S(2,2) 
AG(3,2)=-S(1,1>^S(2,3)+S(1,3)*S(2,1) 
AG(5,5)=S(1,1^S(2,2)-S(  1,2)^S(  1,2) 

DO  1  1=1,3 
DO  1  J=1,5 
AI(I,J)=AG(J,I)/D 

1  CONTINUE 
RETURN 
HID 

SU3I0UT INE  MULT  ( A , B ,  C) 

THIS  ROUTINE  MULTIPLIES  WO  MATRICES 

DIFiEISION  A(3,3),B(  3),C(3) 

DO  1  1=1,3 
SUM=0.0 
DC  2  1  ^ 

2  SUM=SIII+A(I,K)->;B(I0 
C(I)=SUI'l 

1  CONTINUE 
RETURN 
END 

SUBROUTINE  SII'Q(  A,B,N,KS ) 

SOLUTION  OP  THE  LINEAR  ALGEBRAIC  SISTEH  OP  EQUATIONS 

OP  OHE  PORK 

AX=3 

A  -  NXN  mTRIX  OP  COEPPICIENTS 

B  -  VECTOR  POR  THE  RIGKTMD-SIDE  OP  THE  SYSTEM  OP 
EQUATIONS 


o  o  o  o 


202 


N-  NUFiBER  0?  EQUATIONS  AND  VARIABIES 
KS  -  OUTPUT  DIGIT 

0  ?0R  A  NORMAL  SOLUTION 
1  FOR  A  SINGULAR  SET  OF  EQUATIONS 
DIMENSION  A(1),B(1) 

T0L=0.0 
KS=0 
JJ=  -N 

DO  65  J=1,N 

jy=j+ 1 

JJ=JJ+N  +  1 
BIGA=0.0 
IT=JJ-J 
DO  50  I=J,N 
IJ=ir+I 

IF(ABS(BIGA)-ABS(  A(IJ) ) )  20,50,50 
20  BIGA=A(U) 

IM/iX=  I 
30  CONTINUE 

IF(ABS(BIGA)-TOL)  '55,55  ,40 
55  IS=1 
RETUr^J 

AO  I1=J+N''KJ-2) 

IT=MAX-J 
DO  50  K=J,N 
11=1 1+N 
12=11 +IT 
SAVB=A(I1) 

A(I1)=A(I2) 

A(I2)=SA^vE 
50  A(Il)=A(ll)/ErGA 
SAT£=B(mAX) 

B(IMAX)=B(J) 

B(  J)=SAVE/BIGA 
IF  (J-N)  55  ,  70,55 
55  IQ5=N  >' (J-1  ) 

DO  65  K=JY,N 
nJ=IQS+IX 
IT=J-IX 
DO  60  JX=JY,N 
IXJX=1'MKJX-1  )+IX 
JJX=IX  JX-f  IT 

60  A(  I X JX )  =A ( IXJX ) -(A  (  i:gj  )*a(  jjx ) ) 

65  B(IX)=B(n)-(B(J)>A(IXJ)) 

70  NY=N-1 
IT=LWN 

DO  80  J=1,NY 
IA= IT-J 
IB=II-J 
IC=N 

DO  80  K=1,J 
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B(IB)=B(IB)-A(IA)^^(IC) 

IA=IA-N 
80  IC=IC-1 
RETURN 
2ND 

SUBROUTINE  FORM(N,M) 

THIS  ROUTEIE  FORI'IS  THE  COEPHCEINT  MATRICES  USED  TO  SOLVE 
FOR  ELASTIC  BITERIAraNAR  STRESSES 

DIMENSION  THETA(IO) 

COMIi)N/XrjMIT/DC(  2,51 )  .DB(6, 51  )  ,EP 
COMMDN/ROT/Z9 ( 10000 )  ,D(  100  )  ,NA.LF 
C0M^DN/T0F/E5,  T2 ,  V 1 , 3'IX  ,  3'IY  , G ,  CONS  ,  Q ,  A1 ,  SCON 
COI#DN/BO T  /E4 ,  T4  ,  V4  ,  sex  ,  SCY  ,  GB  ,  CONB  ,  Q 1 ,  GC  ,  Ep 
COMON/ADHE  S/TAD ,  GAD 
COIt-DN/BOND/F , P1 , P2, XKF  ,XKUNS  ,XKSTIP 

COMMON /B0ND2 /NE  ,NL  ,1OT  ,XD(  100  ) ,  YC(  100  )  ,XA(  100  )  ,YA(  100  ) , 
INOP(IOO) 

COMMDN/CTOL/TOL  ,NC  ,im  ,TX  ,TY  ,NBC(  100  )  ,LBC 
COMPLEX  Z,CI,P1,G2 
EXTERNAL  F1,G2 
CI=CMPLX(0.  ,1.0) 

A=A1  ! 

0=M/2 

?=N/2 

L=P 

MX=NC*NTl 

K=0 

ASSEI4BLE  RIGHTHAND  SIDE  VECTOR 
PRINT  57 

57  FORMAT(/*  NODE  NUMBER*8X*X-00R*8X*Y-C0R* /) 

PRII7P  56,  (I  ,DC(  1 ,1)  ,DC(  2,1)  ,1=  1,MX) 

56  FORMAT(I10,2E15.3) 

DO  101 1=1, K 
X=DC(1,I) 

Y=DC(2,I) 

Z=X+CI*Y 

CALL  REMro(DU,DV,Z) 

CALL  RErmJ(BU,BV,Z) 

D(l)=nJ-BU 

D(K+I)=DV-F/ 

DO  102  j=i,m:< 

CALIXINTG(Z,J,A,C1,C2,C3,C4,F1,1  ) 

GALL  nNTG(Z,J,A,B1,B2,B3,B4,E1,2  ) 

802  CONTIJJUE 
11=1 
J1=J 

IV1=(J1-1  ><-N+I1 
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Z9(IV1)=C1+B1 

IF(  I.H3.J)Z9(IV1)=Z9(IV1)+TAIiMD 

12=1 

J2=J+L 

IY2=(J2-1  y  11+12 

Z9(IV2)=C24B2 

I3=I+K 

J3=J 

IV3=(J>1  y  11+13 
Z9(IV3)=C3+B3 
14=1 +K 
J4=J+L 

iy4=(J4^1  )*N+I4 

Z9(IV4)=C4+B4 

IF(I4.BQ  .J4)Z9(IV4)=Z9(IV4)+TAI)MD 
102  CONTINUE 
101  CONTINUE 
104  OONTimE 
20  CON  T  DUE 
EETURN 
END 

SUBROUTINE  REMTU(U,V,Z) 

CALCULATE  DISPLACEI-IENTS  IN  THE  TOP  SHEET  DUE  TO  A  RET-IOTE 
STRESS 

OOMMCN /T0P/E3,  T2 ,  VI ,  S'K ,  OT ,  G , OONS  , Q ,  A1 ,  SCON 
COMPLEX  Cl , Z ,  ZI ,  ZB  , CXSR  ,D , DPHI  ,iHI , DOMLBGA 
CI=CMPLX(0.  ,1.) 

ZB=OONJG(Z) 

D=CXSR(Z,A1) 

DPHI=  .  5  *SMY*D- .  25*Z*  ( SMY-affiC ) 

PH  1= .  5  *SH  Y*  Z/D- .  25  *(  SMY-SMX  ) 

D0MEQA= .  5*SITY‘CXSR(  ZB,A1 )+  .25  *(SMY-ayiX)*ZB 
D=  ( Q  *DPH  I  -DOME  GA-  (  Z  -ZB )  *  C  ON  JG  ( PHI ) )  /  ( 2 .  *G ) 

U=REAL(D) 

V=Air4AG(D) 

RETURN 

END 

SUBROUTINE  REI''IBU(U,V, Z) 

CALCUIATES  DISPLACDIENTS  HI  THE  BOTTOM  SHEET  DUE  TO  A 
RH40TE  STRESS 

COMHOtl /EOT /E4,  T4 ,  V4 ,  SCX , SCY , GB  , OONB , Q 1 ,  GC  , E5 

COMPLEX  Z 

X=RIAL(Z) 

Y=A]MAG(Z) 

S12=-V4/E5 

U=  ( SCX/E54B  12*SCY  }*X 
V=  (S12-!^SCX+9SY/E4)*Y 
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RETURN 

MD 

COMELEX  FUNCTION  F1(Z,Z0) 
EVAEUATE  INTEGRAL  OF  B(Z,Z0) 


COMIEEX  Z,Z0,ZB,Z0B,Xft.,XB,XC  ,GXSR,CI  ,D,B,R 
COMPLEX  H,H1,H2.BI 
DIMINSION  N(300) 

C0IWDN/T0P/E3,  T2,  V1 ,  Sa  ,  ffilY , G ,  CONS  ,Q ,  A1 ,  SCON 
COMMON/KK/KCOUNT 

COUNTER  LIMITS  INTEGRATION  INTERVAL  TO  50  STEPS 
USE  comai  to  zero  counter 

XA(Z,Z0)=CL0G((Z*Z0-A**2~GXSR(Z,A)*CXSR(Z0,A)y(Z*Z0+A 
2*^  2- 

2CXSR(  Z,A)*CXSR(  ZO,A)  )*  (-1 .  )) 

XB(  Z , ZO  )  =CL0G (  ( Z*Z0-A*-* 2- CXSR (  Z , A)* CXSR ( ZO , A) ) * (  Z+  ZO  )* 
1*2/ ( 

2 (Z* ZOf A* * 2-  GXS R ( Z, A) * CXSR(  ZO , A) ) *  (  ZO-Z) ** 2 ) ) 

XC ( Z ,  ZO )  =  ( Z-ZO* CXSR(  Z , A)  /CXSR ( ZO ,  A) ) *  ( ZOB-ZO )  /  ( Z**2- 
1ZO**2) 

G1=G 

ZB=OONJG(Z) 

Z0B=C0NJG(Z0) 

CI  =  CMPLX(0.  ,1. ) 

SAVE  PREVIOUS  VALUES  OF  XB  MD  CHECK  PATH  FOR  BRANCH  CUT 
IF  BRANCH  IS  CROSSED  ADD  PROPER  IMAGINARY  TERM  TO  KEEP 
DISFLACEI'IEHTS  CONTINUOUS 

KC0UITT=KC0UNT+1 

K=KOOUNT 

A=A1 

BI  =  .5*(XB ( Z, ZO ) -Qi^XB ( Z,ZOB )  ) +XD(  Z , ZO ) 
X=REAL(XB(Z,ZOB)) 

Y=AIMAG(XB(Z,ZOB)  ) 

1=1 

IP(X.LT.0.AND.Y.GE.0)I=2 

IF(X.LT.0.AND.Y.Iff.0)l=5 

IfcXB(Z,ZOB) 

IF(K.LE.2)  GO  TO  100 

IF(N(K-2).B3  .2.AND.I.B0.3)D=D-6.283*CI 
IF  ( N ( K-2 )  .B3  . 3  .  AND .  I .  BO  .  2  )D=  D+6 . 283  ^I 
100  CONTINUE 
N(K)=I 
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C 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


H1=-Q-“XC(Z,ZO) 

H=  -Q  -  RSAlX  X A  (  Z ,  ZO ) ) 

H2=  .5  *(  Q  C  OI'UG  ( D ) ) 

R=H+H1+H2 

F 1  =R+ XC(  ZB ,  ZO ) + 2*  ( Z  ZO-ZOB^^  Z  B )  /(ZB*  2-  ZO*  *  2 )+  ( Z  -  ZB )  *  B  ( 
1ZB,ZO) 

01  =  5.75*10*^-6 

F1=F1«OOITS 

D1=REA1(D) 

D2=AIM/iG(D) 

RETURN 

END 

SUBROUTINE  ZINTO(Z,I ,  A,S11  ,S12,  S21 ,  S22  ,F1,  IC) 

THIS  ROUTINE  INTEORATES  BE  SU14  AND  DIFFERENCE  OF  TV/0 
COI'ELEZ  FUNCTIONS,  PI  AND  F2,  USED  'TO  EVALUATE  TIE  GREENS 
FUNCTIONS  FOR  DEPLACEMEI-ITS. 

COI-H-'DN/CTOL/TOL  ,NC  , M ,  TK  ,  TV  ,1EC  (  100  ) ,  IBC 
COiaiON/IOC/KCOUNT 

COMPLEX  P 1 , F2,  Z , ZO  ,CI , A1 ,  B1 ,  ZOB 
COMOIT/jEIKIT/DC  (  2 ,  51 ) ,  IB(  6 , 51  ) ,  TF 
EXT3UIAL  F1 
CI=ai?LX(0.  ,  1.) 

X=REAIj(  Z) 

311=0.0 

312=0.0 

321=0.0 

S22=0.0 

Y=AirAG(Z) 

I?  inde:\=i  then  z  lies  outside  integeitioit  patch 

ANDK1  iOINTS  APE  USED  DI  A  SINGLE  INTEGRATION. 

IP  DITJEX=2  THEN  BEN  Z  LIES  V/ITHIN  BE  INTEGRATION  PATCH 
AND  TWO  INTEGRATIONS  ARE  IIADE  V/ITH  K2  POINTS  IN  EACH 
II-ITEEIATION  . 

ICOUNT  CONTROLS  THE  LOOPS  ON  THE  INTEGRATION.  IF  ICOUNT 
EQUALS  0  CONTINUS  OTHEiMISE  DO  SECOND  INTEGRATION. 

XST/iRT=DB(l,l) 

:EINAL=D3(2,I) 

IF  ( X  .IE  . XS Tx\RT  . OR  .X .  GE  . XFI NAL)  INDE:{=  1 

IF(X.GT  .XSTART  .AND  .X. IT  .I{PINAL)INDEX=  2 

IC0UNT=0 

K1=5 

K1  =  5 

IF(niDEX.IQ.1  )2,5 
2  ,E=XSTART 
XU=:iPINAI; 

GO  TO  4 
5  }E=XSBIRT 
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XU=X-TOL 

icouiiT=icomm  1 
GO  TO  4 

5  XU=X PETAL 
:C[,=X+TOL 
ICOUXT=ICOUIIT+1 

^  IF(XF5M;.m)  .A)XPITrAIP=A-TOL 
11=  K1 

ex=(:e-:{l)/(3.*e) 

M=K+1 

DO  7  L=  1  ,M 
XX=(XU-XL)*(L-1  ) 

XO=XX/I5fXI 

AT=L 

X1=AT/2. 

X2=L/2 

IF(1,]33.1  .0R.L.B3.M)16,1S 
16  XF=1 

GO  TO  17 

18  IP(i\BS(X1-X2).IT..0001  )XP=4. 

IP(ABS(X1-X2).GT..0001  )XP=2. 

17  CONTBTUE 

Z0=XOfCI^-0.0 
YO=AIMAG(ZO) 

Z0B=C01IJG(Z0) 

EC0IJIfT=0,0 

CALI  YIl]TG(Z,ZO,I,A,T11  ,T12,  T21 ,  T22  ,P1 ,  IC) 

201  PORKATCA  Y-STRIP  at  XO=*E10.5*T1,T2,T3,T4=*4E11 .3) 
C11=XP^EX*T11 
C12=XP*DX-T12 
C21=CCF^^DX’«-T21 
C22=XF’fD]?=T22 
S11=S114C11 
S12=S12+C12 
521=321+021 
322=522+0  22 
7  con  THEE 

IF(IC0UHT-1  )6,5,6 

6  COIITEIUE 
RETUPIT 
EIlD 

SUBRDUTIIIS  YIIITG(Z,ZO  ,I,A,S11  ,312  ,  321 , 522  ,F1,  IC) 

THI3  ROUTII'IE  iriTEGRATE3  THE  SUM  AlTD  DIFPEREITCE  OF  TWO 
COMPLEX  FUllCTIOHS,  F1  AlTD  F2,  USED  TO  EVALUATE  THE  GREENS 
FUI-3GTIQNS  FOR  DISPLACEIiENTS. 

COMEjEX  H,A2,B2 

COMT'ION/CPARIVS  1 , 32 ,  C 1 1  ,C  12 ,  C 22  , C21 ,  PI ,  P2,  Q4,  Q2 
COMPLEX  Z1  ,Z2,V/1 , V/2 ,G ,GEB  , 31 , 32, C11  ,C12,  C22,C21,  P1,P2,Q 
14,Q2 
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C 

c 

c 

c 

c 

c 

c 


OOMMON/XLIMIT/IX::(2,51  )  ,IB(  6, 51  )  .iT* 

COMMai /CTOL/TOL  ,HC  , NR , TX  , TT  ,1®C (  lOO , IBC 
COKELEX  P1,P2,Z,ZO,CI,A1.B1,ZOB 
COHMaj/ROT/Z9(lOOOO)  ,D(  100  )  ,NA.LF 
CI=CKPLX(0.  ,1  .) 

XD=RE/Jj(ZO) 

X=REAL(Z) 

311=0.0 

312=0.0 

321=0.0 

322=0.0 

y=A]I'IAG(Z) 

B1=DB(3,I) 

B2=DB(4,l) 

C1=DB(5,I) 

C2=DB(6,I) 

Y1=YYX(X0,B1,B2) 

Y2=YYX(X0,C1,C2) 

IC0U1'IT=0 

K1=3 

K1  =  5 

XL=Y1 

XU=Y2 

11  =  1 

IF(n.BQ  .1)G0  TO  4 

IjOGIC  to  STATEI^rEI'IT  FOUR  IS  FOR  THE  NARROV/  STRIP  THAT 
CONTAINS  OHE  SB'IGULARITY 

THESE  STATEOTTS  CAN  BE  USED  TO  INTEGRATE  THIS  STRIP  V/ITH 
ADDITIONAL  LOGIC.  CONSIDER  X  CONSTANT  OVER  THIS  NARROW 
STRIP 


IF  ( Y. IE  . Y 1 . OR  .  Y .  GE  .  Y2  )INDEX=  1 
IF  ( Y .  GT  .  Y 1 .  AND  .  Y .  IT  .  Y2  )INDEX=  2 
IF(  INDEX  .IQ  .1)2,3 

2  :<L=Y1 
XU=Y2 
GO  TO  4 

3  :ai=Yi 
}aJ=Y-TOL 
IC0UNT=IC0UNT+1 
GO  TO  4 

5  XL=Y+TOL 
:{U=Y2 
K1=5 

IC0U1JT=IC0UNT+1 

4  CONTINUE' 

K=K1 

DX=(XU-ZL)/(3.*K) 

M=K+1 

DO  1  Ip=1,M 
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XX=(XD-XL)*(L-1  ) 

YO=XX/K+XL 

AT=L 

X1=AT/2. 

X2=L/2 

IF(L.B3  .1  .OR.L.IQ.M)  16,18 
16  XF=1 

GO  TO  17 

18  IP(ABS(X1-X2).IT..0001)XP=4. 

IF ( ABS  ( XI- X2) . GT  . . 0001  )XP=  2 . 

17  CONTBIUE 

Z0=XO)-CI->^0 
ZOB=CaMJG(ZO) 

IC=1  CRAGf'IED  I'ETAL  SHEET 
10=2  ORTHOTROHC  SOIID  SHEET 

IF(IC.HE.2)30,31 

30  CONTINUE 
A1=F1(Z,Z0) 

B1=F1(Z,Z0B) 

T11=REAL(A1+B1) 

T12=REAL(CI*(A1-B1)) 

T21=AIMAG(A14B1 ) 
T22=AII4AG(CI*(A1-B1)) 

GO  TO  35 

31  CONTINUE 
Z1=X+S1*Y 
Z2=X+S2*Y 
W1=X0+S1*Y0 
V/2=X0  +  S2^Y0 
A1=G(  Z1  ,vn) 

B1=G(Z2,V/2) 

A2=h(zi  ,vn5 
B2=H(Z2,Vf2) 

T 1 1 =2 .  *REAL(F  1*0  11  *A  1+P2*C  21  *B  1 ) 
T12  =  2.*REAI(P1*C  12*A2+P2*C22*B2) 

T21  =2.*REAL(Q4*0  11*A  1-tQ2*C21  *B1 ) 
T22= 2  .*REAI(Q4*C  12  *A24Q  2*0  22*B2 ) 
33  CONTINUE 

IP(IC.NE.3)60,61 

60  FS=1. 

FR=1 . 

GO  TO  62 

61  CONTINUE 
FS=D(I) 

FR=D(I+NALP) 

62  CONTINUE 
P11=XP*FS*])5P^T11 
P12=XF*FR>«I>X*T12 
P21=XP*FS^DX*T21 
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?22=XP-^FR>^DX*T22 
S11=S11+?11 
312=S12+P12 
S21=S21  +r'21 
S22=S22+?22 
1  COIITIrJUE 

IP(IC0I1IT-1  )6,5,6 
6  CX)1TTINUE 
RETUH'I 
HTD 

COMl'LEX  FUIICTION  CXSR(Z,B) 


THIS  FUNCTION  TAKES  THE  SQUARE  ROOT  OF  Z*^2-B**2 
IT  ONLY  RETURNS  THE  FIRST  ROOT,  THE  SECOND  ROOT  CAi\  Bji 
FOUKD  BY  ADDING  IT  /2  TO  THE  TRIP  ARGUILStlT 


COMPLEX  Z,T1,T2,CI 
CI=CI':PIX(0.  ,1.) 

T1=Z-B 

T2=Z+B 
X1=REAL(T1 ) 

Y1=AIM/iG(Tl) 

X2=REAL(T2) 

Y2=AINL&G(T2) 

A1=ATAN2(Y  1,X1  ) 

A2=ATA1I2(Y2,X2) 

S 1 = ( REAL  ( T 1  )■>=•  *  2+AI HAG(  T  1  >  *  2  . 5 

S2=(  REAL  (  T 2 2+A IMAG  ( T2  )**  2  . 5 

SR=(S1^^S2)-->'.5 

ANG=A1+A2  ^  ,  ,  SN 

CXSR=SR=  (COS(AlJG/2  .0)+CI*SIN(ANG/2  .)) 
RETURN 
ET'ID 

COMPLEX  FUI'ICTION  NK(Z,ZO) 


GREENS  FUNCTION  FOR  STRESS  DITENSITY 

COMMON  /TOP/E3,  T2,  V1 ,  SIX  ,  SlY  ,  G,  CONS ,  Q ,  A1 ,  SCON 
COMPLEX  Z ,  ZO  ,  ZB  ,  a)3  ,  CXS R ,  D ,  F 
A-A1 

ZB=OONJG(Z) 

ZOB=OONJG(ZO) 

D=(  ZO*''ZOB-2  .*ZO''^*^'-2+A-^  -  2)/((ZO^'^'2--A>^‘^^  2  )^CXSR(  ZO,A) ) 
F=0xCIvSR(a)B,A)/(Z0B-x-:'^2-A^*2) 

. 3*(F+D ) / (  6 . 283 *(  1 . 4Q  )  ^-  T2 ) 

RETURN 

j3TD 

FUNCTION  Yy(T) 

C 

C  T  IS  A  DUT<FiY  PARAMETER  FOR  XO 
C  CALCUI..ATE  YO  FOR  A  Gr/EN  XO 
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COmOH/TOP/E3,  T2,  V1 ,  S'lX  ,,34Y  ,G, COWS  ,Q, A1 ,  SCON 
COMON/BOND/F ,  P1 ,  P2,  ^QCF  ,XKUNS  ,XKSTIPF 
IF(T/A1.IuC.O)2, 1 
1  CONTIIiaE 

Yr=I^  (  (  1-  (T /A  1  )■**  P 1 )«  *  ( 1  /P2 ))  *A  1 
GO  TO  3 
2  yy=F<Ai 
3  CONTINUE 
PETUM 
END 

COMPEX  FUITCTION  B  (  Z  ,  ZO ) 

CONMON/TO P/E3 ,  T 2,  V 1 , 34X  ,  3IY  ,  G ,  CONS  ,  Q ,  A1 ,  SCON 
COKPTEX  X1  ,X2 ,X3  ,X4,X5  ,CXSR,ZOB  ,  ZD  ,  Z,X6,  X7,  X8,  X9 
Q1=Q 
A=A1 

ZOB  =  COIIJG(ZO) 

X1=CXSR(Z,A) 

X2=CXSR(ZO,A) 

X3=1./(Z*"'2-ZO->:'5^2) 

X4=1./(Z-^-’'2-Z0B-x-2) 

X5=1./(Z-ZO  )•>«••>:  2 
X6=1./(Z+ZO)**2 

X7=X1*(  (-4*ZO*X3)+2.'^Q1*ZOB*X4+(Z-ZOB)*X5-  (Z+Z0B)*X6) 
X8=(  ZO-ZOB )  >^(  (  A^->-2-  Z'^ZO  X5-  (A^  2+Z^  ZO  )*X6)  /]{2+2. *X  2* 

1X3*Z 

X9=  -2 .  1  *Z^^  C  XSR  (  Z  OB ,  A)  *  X4 

B=  ( X7+X  8+X  9 )  /  ( 2 . 0^-X  1 ) 

RETURN 

END 

SUBROUTINE  CPAR 

CALCULATE  PARAHTEPS  FOR  COMPLEX  VARIABLES  IN  ORTHOTROPIC 
iUTALYSIS. 

COMPEX  Cl  ,31,S1B,S2,S2B,C11  ,C12,  C21 ,  C22  ,P1,P2,  Q4,  Q2,D 
1,U1,U2 
COMILEX  CC 

C0IiI40N/B0T/E4,  T4,  V4,  SCX  ,  SCY  , GB  ,  CONB ,  Q1 ,  GG  ,  E5 
C0I'EI0N/Cl-  ARri/31,S2,C11  ,C12 ,  C22  ,C21 ,  P1 ,  P2,  Q4,  Q2 
CI  =  CELX(0.0, 1.0) 

EX=E5 

IY=E4 

VX=V4 

GXY=GC 

B=EX/G:CY-2*VX 

C=EX/Ey 

CC=CM'LX(C,0  .0) 

VY=7X'->:  EY/E:[ 

D=CSQRT(B-”-2-4.*CC) 

U1=-.5*(B-D) 
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U2=-.5*(B+D) 

S1=CSQRT(U1) 

IF(AIl'![AG(Sl).nD.O  )S1=-S1 
S2=CSQRT(U2) 

IP(ABIAG(S2).ni:  .0)S2=-S2 
S1B=C0NJG(S1) 

S2B=C01UG(S2) 

IP(B^^2-4.*C)  1,2,3 

1  CONTINUE 
PRINT  4 

4  FORFiAT(*  ROOTS  ARE  COMPLEX  U1=A+IB,  U2=-A^IB*/) 

GO  TO  5 

2  CON  T  HUE 

PRINT  6  „ 

6  PO?J-iAT(*  lOOTS  .^RE  PURE  MAG  I  NARY  AND  EQUAL*/) 

GO  T05 

5  B1=AIMAG(S1) 

D1=A.IKAG(S2) 

P1=(S1*->^2-VX)/EX 
P2=(S2**  2-'VX)/EX 
Q2=(-VY*S2+1./S2)/IY 
Q4=(-Vy*S1+1./Sl)/E5r 

ICO  PORMAT(4E12.3) 

COMP1= 1  .  /( 12 . 566-’^  4*(B  1**'  2-D1**  2 )) 

C 11  =(D  1**  2*VY+ 1  )*B  1*COMP1 
C 12  =(  VX+D 1**  2  )*COKP1 
C21  =-D1*(  VY*B1**  2+ 1  )*C0MP1 
C  22= - (B 1  **  2+VX )  *C  OMP1 
GO  TO  7 
5  IRINT  8 

8  P0RI'1AT(/*  ERBDR-ERROR  ORTHOTROPIC  ANALYSIS  IS  NOT* 

1*  DEFINED  */) 

7  CONTINUE 
RETUHT 

mD 

SUBROUTINE  GRID 

GENERATE  MESH  USED  TO  DISCRETIZE  INTERLAMINAR  STRESSES. 

COMMON /BOND2An^E,NL  ,NT  ,XC(  100  )  ,YC(  100  )  ,XA(  100  )  ,YA(  100  ) , 
INOP  (100) 

COMMON  /CT 0 L/TOL ,  NC  ,  NR ,  TX ,  TY ,  NBC  ( 100  ) ,  LBC 
C01M0IT/T0P/E3,  T2,  V 1 , 3'IX , S-IY , G ,  OONS , Q ,  A1 ,  SCON 
COMMQI  /:OLIM  IT/DC  (  2 ,  51  ) ,  IB ( 6 , 51  ) ,  LP 
COMMON/BOND/P ,  PI ,  P2,  XKP  ,XKUNS  ,XICSTIF 
P=P1 

TXX=2.*T0L 

DX=TX/NC 

DY=TY/NR 

LBC=0 

MX=NC*NR 


oooo  ooooo 
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DO  11  1=1, MX 
11  KBC(I)=0 

DO  100  J=1,NR 
DO  100  1=1, NO 
MT=iro*  (J-1  )+I 
AG=A1+DY*(J-1  )+DY/2. 

YI  =  P*A1 

BC=YI+DY*(  J-1  )+DY/2  . 

DO  ( 1  ,MT  )=DX/2  .+  ( 1-1  )*DX 

x=dcO,mt) 

Y=YYD(X,BC,AC  ,P) 

DC(2,MT)=Y 
DB(1,MT)=DX:<(I-1 ) 

DB(2,KT)=DX*I 
DB(3,ro)=I>Y*(J-1  )+YI 
DB(4,MT)=DY^(  J-1  )+A1 
DB(5,I«)=Dr*J+YI 
DB(6,ME)=IY*J+A1 
100  COKTIIIUE 
110  K)RMAT(I10,2E15.3) 

120  F0RMA[P(I10,6E12 .3) 

RETURN 

EfJD 

SUBRDUTDTE  XENC  (SUI>I  ,XNCY  ,mS  ,DPS ) 


THIS  ROUTINE  INCREMENTS  THE  CRACK  LENGTH  MD  DEBOND  SHAPE 
FOR  NCY  CYLES 

COM-DN/XLIMIT  /DC  (  2 , 51  ) ,  IB(  6 , 51  )  ,  EP 
CCMI^N/CTOL/TOL  ,NC  ,HR  ,TX  ,TI  ,NBC(  100  ) ,  LBC 
COMMON/BOND/P ,  ?1 ,  P2 ,  XKF  ,XKUNS  ,XKSTIFF 
COMMON /TO I7E3,  T2,  VI ,  S^IX  ,a:'TY  ,G,OONS  ,Q,  A1 ,  SCON 
COMMaJ  /DD T  /E4 ,  T4 ,  V4 ,  sex  ,  SCY  ,  GB  ,  CCNB , Q 1 ,  GC  ,  E5 
COmON/XY/XD  ,YD 

DIMZtTSION  TiR  20  )  ,  XN(  20  ) , BT  (  20 , 2 ) ,  Cl  (  20  ) , D(  20  )  ,CR(  2 , 2 ) 

B1=A1*F 

R=.01 

DA = 3 . 2 2E-1 4  *XKS  T IFP**3 . 38 

CFAC=1.79E-14 

CFA0=3.36E-14 

DA =CF  AC^  XKST I FF^*  3 . 38 

DA=DA/((  1.-R)  *56000. -XKSTIFF) 

DF=3. 158E-05*SUM**3. 616 
DAS=D& 

DFS=DF 


DETERMDTE  HOW  MANY  CYCLES  REQUIRED  FOR  EITHER  A  CRACK  OR 
DEBOND  EXTENSION  OP  .1  INCHES.  THEN  USE  SMALLEST  VALUE 
AS  THE  INCREMENT  OF  APPLIED  LOAD  CYCLES 


XNCRACK=.10/riA 


oooo  oooo  ooo 
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jaJB01®= .  10/IF 
XNB0ND=.20/IF 
X1JCY=XNCRACK 

IP  ( XNGRA CE .  GT  .  XNBO ND  )XN CY=XNB0ND 
rA=XIICY*DA 
A2=A1+I»A 
DF=DF*XIJCY 
B2=B1+DF 
?1  =  2 
A1=A2 
F=B2/A2 
RBTURIT 

END  _ _ 

SUBHOUTBIE  VErO:  ( Z  ,SXX ,  SXY  ,  SXXY  ,SYX  ,  SYT  ,  SYXY  ,MX  ,KTYPJi ) 

BITEGRATE  (FEENS  FUNCTIONS  FOR  STRESSES 

COMMON/CTO L/TOL,NC  ,NR  ,TX  ,TY  ,NBC(  100  )  ,IBC 
COMPLEX  Z 

z=  z 

MX=MX 

J=J 

sxx=o. 

SXY=0. 

SXXY=0. 

SYX=0. 

SYXY=0. 

SIY=0.0 
DO  1  J=1,I«IX 

CALL  VXmTG  (J  ,Z,S1,S2,S3,S4,S5,S6,r'IX,KTYPi) 

FIRST  INDEX  INDICATES  LOAD  DIRECTION,  SECOND  STRESS 
DIRECTION 

3XX=S1+SXX 

SXY=S5+SXY 

SXXY=S5+SXXY 

SYX=S2+SYX 

SrY=S4+SYY 

sy:{y=s6+syxy 

1  CONTINUE 
RETURN 

END  V 

SUBROUTDTE  VX I N TG (I ,  Z ,  S 1 ,  S2 ,  S3 ,  S4 ,  S5 ,  S6 ,  M ,  KT YPE ) 

DETERl'IINE  STRESSES  M  ADHERENDS  DUE  TO  INTERLMINAR 
STRESSES. 

COMOK/XIIMIT  /DC  (  2 , 51 ) ;  IB ( 6 ,  51  ) ,  EP 

COMI>DN/RO  T /Z9  ( 1000 0  )  ,D(  100  )  , NALF 

COMMCN  /TO  F/E3 ,  T2 ,  V 1 , 34X  ,  3'!Y  ,  G ,  CONS , Q ,  A 1 ,  SCON 
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C 

C 

C 

C 

C 

C 

C 


COMIIEX  Z  ,  ZB ,  20  ,  Z03 , B ,  IBZ ,  B 1 ,  DB 1 B ,  31 B ,  HB 1 ,  GS ,  G6 ,  G7 ,  G8 

COMHOK  /CTC IVTOL , IIC  ,  IIR  ,  K  ,  TI , IIBC  ( 100  , LBC 

COMPLEX  Cl 

C0K1B3X  P 

1-^Z 

/:=REAL(  Z) 

31=0 

32=0 

35=0 

34=0 

35=0 

36=0 

ZB=OOITJG(Z) 


IP  IITIiI;X=1  THEN  2  IIS  OUTSIDE  THE  INTEGRATIGIT  PATCH 
AIID  HI  JOINTS  ARE  USED  IN  A  SINGLE  INTEGRATION 
I?  IjDSX=2  TIIEN  Z  LIES  ’\TTHIN  THE  INTEGRATION  PATCH  AND 
Ti.7C  INTEGRATIONS  ARE  I-IEADE  lACH  VIITH  K2  INTEGRATION  POINT 
IF  ICGUNT=0  THEN  CONTINUE  OTHERUISS  SECOND  INTSGilATION 


>START  =  DB(1 ,1) 

X?INAL=DD(2,I) 

IF ( X .IE  . XS TART  . OR  . X . GE  . XPI NAL)  INDEX=  1 

IF(X.GT  .IGTART  .AND.X.IT  .XFINAL)INDEX=2 

K1=3 

K1=5 

IC0UNT=0 

IF(INDEX.EQ  .1)2,  5 

2  ZL=X3TART 
J'IJ=XPINAL 
GO  TO  4 

3  XI.=XSTART 
:-H=X-TC'L 
IC0UNT=IC0UNT+1 
GO  TO  4 

5  :<U=X?IiUvL 
/jj =icf  T  0  j-i 
IC0UI;]T=ICC'UNT+1 

4  CONTINUE 

IFd'EINAL.EQ  .A1  )XFINAL=A1- TOL 
K=K1 

DX=(xu-:c)/(3.-='K) 

M=K+1 

ci=an-LX(o.  ,1.0) 

DC  7  L=1,M 
:c{=(xu-xr.)^-(L-i ) 
xo=:g{/k+xj., 

■ro  =  YY(X0) 

ZO  =  XO+CFYO 
ZOB=CCITJG(ZO) 

AT=I. 
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X1=AT/2. 

X2=L/2 

IF(L.}i)3.1  .OR.L.HJ.M)  16,18 

16  ZI’=1 

GO  TO  17 

18  I? ( /iBS  (X 1- X2 ) .  liT  . .  0 0 01  )XF=  4 
IP(ABS(X1-  X2)  .GT..0001  )XP=2 

17  CX)NTIi]IiE 

CALL  7YI IT TG(  Z , ZO  ,  I ,  TXX  ,  E^Y ,  TXX  ,  TYY  ,  TUX ,  TXYY  ,KIYPS  ) 

S1=XP^E{->:t:{X+S1 

S2=XF->'IIv-TXY+S2 

33=X?*IK*E£X+S? 

34=X?->'-Il>:*'rYY+34 

S5=xp*ix-t:yx+s‘> 

SG=XP*IK-xTXYY+S6 
7  COITTII'TUS 

IF(IC0U1OT-1  )6,5,6 
6  COIITINUE 
RETUITIT 
Si'TD 

SUBiTOUTLTE  VYIITTG ( Z,ZO  ,  I , S7,  38,  S5,  34,  35,  36,  KTYPE ) 

C 

C  INTEGRATE  GREERS  FUNCTIONS  FOR  STRESSES 
C 

COMMON /Cl’ARI'l/S  1 , 32,  C 11  ,C  12 ,  C22 .  C21 ,  PI ,  P2,  Q4,  Q2 
COHMaT /XL  IMIT  /EC  (  2 , 51  )  ,  NB  (  6 , 51  ) ,  PF 

COMPLEX  Z1,Z2,W1,W2,  GB,S1,  S2,C11  ,C12,  C22  ,C21,  PI  ,P2,  Q 
14, Q2 

COr&DN/TCP/EJ,  T2,  VI ,  3>IX  , , G ,  CONS  ,Q , A1 ,  SCON 
C0I'F'DR/R0T/Z9 ( 10000  )  ,D(  100  )  ,NALF 

COMPLEX  Z,Z3,Z0,Z0B,B,IBZ,B1,nB1B,B1B,rB1,G5,G6,G7,G8 
COKPLSG  CI,G1A 

COMPLSTC  G  1 ,  G2 ,  G3 ,  G4 ,  H 1 ,  H2 ,  W 1 B ,  N2B 
COMMCSI  /CT 0  L/TO  L ,  NC  ,  NR ,  TK  ,  TiT ,  NB C  ( 100  ) ,  IBC 
CCMILEX  P 
X=REiff,(  Z) 

Y=AIMAG(Z) 

r=z 

Q1=Q  . 

ZB=Ca'TJG(Z) 

37=0. 

38=0.0 

55=0 

34=0 

35=0 

36=0 

rE{=ITR*NC  ■ 


DB  Giv: 


?ARAI-i 


)TERS 


!R  AND  IDV/SR  GRID  BOUNDARIES. 


B1=DB(5,I) 
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B2=IIB(4,I) 

C1=DBf  5,1) 

C2=DB(6,I) 

XO=REAL(ZO) 

Y2=YYX(X0,C1,C2) 

Y1=YYX(:®  ,B1,B2) 
p:i=5 
K1  =  3 
.'CL=Y1 
:<n=Y2 
K=K1 
C 

C  SEE  YIIITG  FOR  LOGIC  FOR  INTEGRATION  OF  NARROV/  STRIP  THAT 
C  CONTADIS  THE  SEIGULARITY.  THAT  STRIP  IS  NEGLECTED  HERE 
C 

DX=(XU-XL)/(3.-xK) 

M=K+1 

CI=CMPLX(0.  ,1.0  ) 

DO  7  L=1,M 
XX=(XU-XL)*(L-1  ) 

YO=XX/K+XL 

ZO=XO+CI^YO 

ZOB=CONJG(ZO) 

AT=L 

X1^T/2. 

X2=L/2 

IF(L.E3  .1  .0R.L.B3  .M)  16, 18 

16  :<F=1 

GO  TO  17 

18  IF(ABS(X1-X2).LT..0001  )XF=4 
IF(ABSU1-X2).GT..0001  )XP=2 

17  OONTIWUE 
IPCKTXPE.SO.I  )30,31 

30  CONTINUE 
B1=B(Z,Z0) 

B1B=B(  Z,ZOB) 

DB1=DBZ(S,Z0) 

DB1B=DBZ(Z,Z0B) 

G1=-4*(Z0/(Z^-*2-a)**2)+Z0B/(Z-:<*2-Z0B**2)) 

G1=REAL(G1) 

G 1 7b  -4  *(  ZO/  (Z*  *2-  ZO*  *2 )  -  ZOB /  (Z  ^  * 2-  ZOB** 2  ) ) 
G1A=REAL(CI^G1A) 

G2=2->^  1*Z0B/(Z*^-2-  Z0B**2) 

G2=G2-  ((ZOBfZB)/(Z+ZO)**2+(ZOB-ZB)/(Z-ZO)**2) 
G5=2*Q1*ZO/(Z--x2-ZO**2) 

G3=G  3-  ( (ZO+ZB )  / (Z+  ZOB )  **2+(  ZO-ZB )  /  (Z-ZOB ) 2) 
G5=B1+B1B 
G6=CI*(B1-B1B) 

67=(ZB-Z)->'(DB1+DB1B) 

G8=0 1  *  (Z  B-Z)  -  ( DB 1-  DB 1 B) 

T1S=RE7\]7(G1-  G2-G3) 
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T2S=REAL(G1A-CI*(G2-G3)) 

T3S=REAL(GHG24G3) 

T4S=REAL(G1A+CI*(G2-  G3)) 
T5S=AiriAG(G1+G2-tG3) 

T63=  A IMAG  ( G 1 A+ C I  *  ( G2- G3 ) ) 

T 1 1  =  ( -RKAL  ( 2*G  5- G7  HT  1 S)  t  SCON 
T22 = ( -REAL  ( 2*G  6-  G8)+T  2S)  *  S  CON 
T33=  ( -REAL  ( 2*G  5-^0  7  )+T  3S)  *  SCON 
T44=( -REAL  ( 2*G  6+G  8  )+T  4  S)  *  SCON 
T55=(AimG(-G7)+T5S)*SCON 
T66 = ( AllAG  ( -G8  K  T  6  S)  S  CON 

CHAITGE  SIGNS  TO  ACCOUNT  FORNEGATIVE  BODY  FORCES 


T11=-T11 
T22=-T22 
T53=-T33 
T44=-T44 
T55=-T55 
T66=-T66 
GO  TO  32 
31  CONTINUE 
Z1=X+S1*Y 
Z2=X+S2*Y 
V/1=XO+S1^^YO 


v/2=XO+S2-^YO 

V/1B=C0I']JG(W1) 

W2B=C0NJG(VA2)  ,  .. 

G 1  =  2*0.'/ 1  /  (  Z 1  2-  Vn  **  2 )+ V/ 1 B/  ( Z 1  **  2-  ¥1  ^  *  2 ) ) 

G2=  2*0/2/ (Z 2**^  2-  V/2** 2 )+¥2B/  (Z2*x  2-W2B** 2 )) 
H1  =  2*(W1/(Z  1**  2-  V/1**-  2)-V/1B/  (zi*->‘  2-V/1B^*2)) 
H2=  2*(  W  2/  ( Z  2*^  2-  V/2**  2 )  -V/2B/  ( Z  2**  2-  ¥2B>^- *  2 ) ) 

T 11  =  2*RS AL(  S1**  2-'*C  1 1  *G  1+S 2**  2*0  21  ^2 ) 

T22  =  2*REAL(CI*  (SI**  2*C  12*H1+S2*><-  2*0  22*H2 )) 
T33=2*REAL(C11*G1+C21*G2) 

T44  =  2*RE  A  L(  Cl*  (  C 12  *H  1  ■+€  22  *H  2  )  ) 
T55=-2*IUCAL(S1*C11*G1+S2*C21*G2) 

T66  =  -2  *RE  A L  (C I  *  ( S  1*C  1 2  *H  1 +S  2 *C  22  *n  2 ) ) 


CAICULATE  STRESSES  IN  ORTHO TROLIC  SHEET 


32  CONTINUE 

37  =XF  *  DX*  T 1 1  *D  ( I ) + S7 
S8=XF*DX*  T22*D(  I+MX)+S8 
S3=XF* DX* T53*D ( I )+S3 
S  4  =XF  *  DX^-  T  44*D  ( I  )  +  S4 

55  =XF*  DX*  T55  *D ( I )  +S  5 

56  =XF  *  DX*  T66  *D  ( I  +I‘'K )  +  S6 
7  OONTIinJE 

RETURf'I 

HID 
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COMIIiSX  FUNCTION  DBZ(Z,ZO) 

COMO  N/TO  P/E3 ,  T2 ,  V 1 ,  SMX ,  SM  Y ,  G ,  CO  NS ,  Q ,  A 1 ,  SC  ON 

THIS  FUNCTION  IS  THE  DERIVATIVE  OF  B(Z,ZO)  WITH  RESPECT 
TO  Z 

COICELEX  ZB,CXSR 
COMPLEX  Z,ZO,B1  ,B2,B3,ZOB 
COMITiEX  XIZO  ,XI  ZOB , XI  Z 
ZB=CCN  JG(Z) 

Z03=C0NJG(ZC) 

A=A1 

XIZ=CXSR(Z,A) 

XI  ZO=CXSR(ZO,A) 

XI  Z0B=O:SR(Z0B,A) 

B 1  =  (  -4  *Q  -5^  Z*  ZOB  /  ( Z  *  *2-  ZO  B**  2 )  *  *  2-  4  *Z^  ( 3*  Z  0  2*  Z  OBf  Z*  *  2 

1*ZOB-4*ZO^-^  3 

2)/((ZO-Z)'<->^3*(  ZO+Z)**3))/2 

B2=- ( ( Q*X I ZOB-J^  (A*>s  2*Z OB^-  2-2  *Z^'  -^4+A •!<*2*Z * *2 )  / ( ZOB>^-^ 4-2 
1*Z-x*2*Z03 

2**  2+Z«^4  )-XI  ZQ:'  ( A^'^'  2*ZO^  2-2*Z*  H+A  '*2*Z>^*2)/(ZO-x*4-2->' 
1Z.x*2*ZO’'^2 

2+Z'“*4))  /(Z^=--2-A^  *2)+  ( (Z0**2+Z*Z0-2*A^  x2)/(20-Z)^-  ^3+(  ZO 
1-X  X2-Z-X  ZO 

>2*A*^2)/(SO+Z  )^*3+Z^  (( A^-^  2-  Z^  SO)/(ZO-Z)*->=  2-  (Z>  ZO+A->^-*2 
1)/(Z0+Z)**2 

4  2- A-  * 2 ) )  (ZO-ZO  B)/(  2*X I  ZO  )  )/XI  Z 

DBZ=B14B2 

RETURN 
EICD 

SUBRC^  UTI NE  RSSIGM  (Z,S11,S22,S12,  KT  YPS  ,  ST  RESSM ,  STRESSC  ) 

C 

C  COr-TUTS  STRE'SiiES  IN  ADIIERENDS  DUE  TO'  RiI>iOTE  STRESSES 

C 

C0I310N/T0I7E3  ,T2  ,V1  ,SKX,SI'i  Y,G,CONS,Q  ,A1  ,SCOi: 

DI  Mrll  S 10 1:  S  T  RES  SI--(  3  )  ,  ST  RES  S  C  ( 3  ) 

CCMIIEX  Z,CXS.R,I'HI  ,0MEGA3,D1'HI,XK;  ,ZB  ,SI 
IF(l'irYPE.EQ.l)l  ,2 
1  CONTINUE 
ZB=CGNJG(Z) 

IHI=ST?ESJf  i(  2)”'  (Z/CXSR(  Z  ,A1 ))  /2 ST;ISS3M(  2)-STRES 

isr:(i)) 

0MEGAB=STRESSK(2)-><  (Z3/CXSR(  ZB,A1 ))  /2  .+  .25*(  STRESSMCZ)- 
1  ST  IISS  3*^  ( 1 ) ) 

DI  KI=  -A  1  - 2 / ( C  XSR  (  Z ,  A 1  )^  (  2-  A 1  2 ) ) 

Di-  H 1=  DIH  1*3  T  RE  S  S'-;  ( 2  )  /2  . 

SI  =  CO  IT JG  ( aiE  GAB  )  -PII  I- Z*  DHil 
^a'HI=2*R£AL(iHl) 

XK=XlHn+  ZB*DPIII+  SI 
S 1 1  =REAI,  ( ICLHI-  ( 2B*DPHI+  SI)) 

S22=REAI,  ( XIC ) 
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si2=aim;ig(xi:) 

GC)  TO 

2  311=3TRESSC(1) 

322=STRESSG(2) 

S12=STRSS3C(3) 

3  CONOIENUS 
RETIRH 
EIJD 

COKi^LEX  PUIfCTION  G(Z,V/) 

COMPLEX  Z,ZB,V/,v'ffl 
ZB  =  OOKJG(Z) 

Vffi=OOMJG(V/) 

G=  C  LO  G  ( ;/-Z )  -CL  OG(  Z+¥B  )  -CLO G ( Z+  ¥) +CLOG  ( Vffi-Z) 

RSTUIGT 

STD 

CCiaLEX  FUlvCTION  H(Z,¥) 

COMPLEX  Z,ZB,¥,VB  ,CI 

Cl  =a'IPIX  (0.0,1  .0) 

ZB=00IJJG(Z) 

V/B=00KJG(V/) 

H=  CLOG(  V/- Z)  +  CLO  G ( Z+  VS  ) -CLO G( Z+ V/) -CLO G (WB-Z) 

H=CI>>H 

RETURN 

END 

SUBRDUTBTE  PDA  STIC  ( SYIELD  ,GAD2  ,  S  ,DR  ) 

C 

C  PSRFCRI-:  INCRJDISNTAL  PLASTIC  ANAI^YSIS 

C 

C0MM0N/R0T/Z9(  10000)  ,D(  100  )  ,RALF 
COKMOIT/AD  HES /TA  D ,  GA  D 

DIMSiSION  F(  100  )  ,G(  100  )  ,IFt(  100  )  ,KOI'(  100  )  ,I'IYIEID(  100  ) 

ac=o 

c 

C  COI'LUTE  HELD  STEiESS  FOR  EAST  POINT 
C  CHOOSE  CRITICAL  ELEKaiT 

c 

K=2*ITAIiP 

NSQ=N^':: 

SUM=0 

DC  55  1=1, NAIF 
F(I)=0.0 
F(l+NAiF)=0.0 
55  N0P(I)=0 

FC I  =993999  999  9  939. 

12  FCC=99  993  9  9  99. 

K}C=iZ+1 
DC  11=1,  HALF 

IF(N0P(I).SQ  .1)G0  TO  1 
C 

c 


ADD  LOGIC  TO  SKIP  ALL  YIELDED  EIEKEIIS 


ooo  ooo  ooo 
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A=D(I  )**2+D(  I+MLF)**2 
3=  2* ( F ( I ) D (I) + ? ( I+FAIi? ) ( 1+ NAI.F ) ) 

C=F(I)**2+F(  I+MLF)  >(*2-SY'IELD*-‘2/5. 

TE S T  =  ( -Bf  SQ  RT(  B** 2-4  *A* C  ) )  /  ( 2 .  *A ) 
IF(TE3T.LT.FCC)2,3 

2  K=I 
FCC=TE3r 

3  CONTIIJUE 
1  CONTINUE 

IF(KK.B2  .1)FCI=FCC 
LL=K-1 

SUM1=SUM 
EUI-1=SUE+FCC 
IF  (  SUE. GT.  3)5,6 

6  IF(FCC.B3  .9  9  9999  9  9  9. )G0  TO  18 

N0P(K)  =  1 
NYIELD(I30=F 

SAVE  STRESS  IN  EAQI  ELEEENT  ATYIELDniG 

DC  8  1=1  ,N 
F(I)=PCC*D(I)+F(I) 

G(  I  )=FCC  -D(  I ) 

8  CONTINUE 

MODIFY  EQUATION  SET  IC'R  YIELDING  OF  CRITICAL  ELEI'IENT 

I1=(K-1  )*N+K 
I4=K+NAI.F 
I2=(I4-1  )»N+I4 

Z9(1 1  )=  Z9(1 1  )+ TAD>^  ( 1 .  /GAD2-1 . /GAD) 

Z9 (1 2 )=  Z9  (1 2 )+TAD^  (  1 .  /GAD2- 1 .  /GAD ) 

USING  AN  ITERATIVE  HETNOD  UE'DATE  ELEKE^IT  STRESSES 

GALL  GAUSS  (DR) 

IF(KK.EQ  .J:mLF)G0  TO  18 
GO  TO  12 

5  IF(H:.E2  .1  )15 , 14 
15  IRINT  15 
3UM=0. 

15  F0R1«1AT(/*  THE  SOLUTION  IS  COMPLETELY  ELA.STIC*/) 
GO  TO  16 

14  CONTINUE 
EE=K-1 
SUM=SUM1 
18  CONTINUE. 

PRINT  31,KK,SUM 

31  K)REAT(/I10*  ELEMENTS  NAVE  YIELDED  AT  *E10.5/) 

16  CONTINUE 

1001  FORMAT (GEIO.  3) 


o  o  o  o 
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DO  26  1=1,  mLF 

D(l)  =  (S-SaM)-D(l)+F(l)  ^  ^ 

D(  I +KALF)  =  ( S-SUM  )*D ( I+NAIiP )  +P(  I+NALF ) 

26  OONTDIUE 
IF(KK.GT.1  )25,29 

25  PRINT  27 

27  K)RKAT(/*  YIHLiDSD  EIEKENTS  */) 

IRIITT  28,(NYIEID(I),I=1,KK) 

28  F0RMAT(12I10  ) 

29  OONTIIIUE 

PRINT  30,FCI  ^ 

50  FORMAT(-''  THE  YIELD  mCROSCOPIC  STRESS  IS  *E11.3/) 

RETUHfJ 

END 

SUBROUTINE  GAUSS (DD) 

GAUSS  SEIDEL  METHOD  FOR  SOLVING  INCREMENTAL 
H-ASTTC  SOIUTiai 

DIMENSION  ID(  100  ) ,  ASAVE(  100  ) 

COMI^N/ADHES/TAD ,  GAD 
COMMON /R0T/Z9(  10000)  ,D(  100  )  ,NALF 
EPS=.005 
N=2*NAIF 
NSQ=1P'-N 
ITMAX=20 

DO  35  1=1  ,N 
K=(I-1)*N+I 
ASTAR=Z9(K) 

ASAVE(I)=ASTAR 
DO  3  J=1,N 
II=(J-1)*N+I 
Z9(II  )=Z9(II  )/ASTAR 
3  CONTUraE 

DD(I)=DD(I)/ASTAR 
33  CONTINUE 

DO  9  ITER=1,ITMAX 
KFLAG=  1 
DO  7  1=1,  N 
:{STAR=D(I  ) 

D(I)=DD(I) 

DO  5  J=1  ,N 
II=(J-1  >><11+1 
IF  (I  .133.  J)  GO  TO  5 
D(I)=D(  I)-Z9(II)*D(  J) 

5  CONTINUE 

IF(ABS(XSTAR-D(l)).IiE.EPS)GO  TO  7 
KPLAG=0 
7  CONTINUE 

IF(KFLAG.NE.1  )G0  TO  9 
GO  TO  1 


o  o  o  o 


9  CONTINUE 
PRINT  204 
1  CONTDIUE 

RECOIBTRUGT  Z9  MD  DD  MATRIX  K)R  USE  ON  EUOUTRE 
ELEI'EIITS  THAT  HELD 


DO  350  1=1,  N 
K=(I-1  )*]A  +  1 
ASTfi.R=ASAVE(l) 

DO  30  J=1  ,N 

n=(J-1  >11+1 

Z9(II)=Z9(II)^ASTAR 
30  OONTINUE 

DD(I)=DD(I)*ASTA.R 
330  CONTINUE 

204  P0RI-'AT(A  CAUTION  THE  GAUSS  SEIDEL  DI 
RETUim 
HID 

FUNCTION  yYX(X,B,A) 

COMMON  /BOND/P ,  PI ,  P  2 ,  XKF  ,  XKUN  S ,  XICS  T  IFF 
IF(X.B3.0)G0  TO  4 
IP(X/A.(2).1  )1,2 

2  YyX=B*^(  (1.-(X/A)^  *P1)*  ’^(1./P2)) 

GO  TO  5 

1  Yyx=o. 

GO  TO  3 

4  YYX=B 

3  CONTINUE 
RETURtT 
END 

FUNCT ION  YYD  (  X , B ,  A  , P) 

IF(X/A.GB.1  )1,2 

2  yYD=B^(  (  1 .  -d/A)* *P  ) *H  1 .  /P  )  ) 

GO  TO  3 

1  YyD=0. 

3  CONTINUE 
RETURN 
EIID 


DID  NOT  COinrSRGE*/) 
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